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CHAPTER  I 
INTRODUCTION 


1.1  INTRODUCTION 

Noise  and  vibration  generated  by  the  rotating  mechanical  equipment  including 
geared  drives  have  always  been  a  problem  in  the  implementation  of  new  technology  in 
automobiles,  rotorcrafts  and  industrial  machines  [1-6].  Recently,  the  need  for  reliable 
vibration/noise  prediction  methods  have  been  found  to  be  crucial  as  faster  and  lighter 
machines  are  being  designed  [6-8].  In  most  of  these  rotating  systems,  structure-borne 
noise  paths  through  bearings,  which  support  the  rotating  shafts  on  flexible  or  rigid 
casings,  are  dominant  [5,6,9].  Hence,  in  order  to  obtain  reliable  mathematical 
prediction  of  the  overall  dynamic  system,  a  complete  understanding  of  the  vibration 
transmission  mechanism  through  bearings,  and  the  role  of  bearings  as  a  dynamic 
coupler  between  the  shaft  and  casing,  is  critical. 

Current  bearing  models,  based  on  ideal  boundary  condition  or  purely  translational 
stiffness  element  description,  cannot  explain  how  the  vibratory  motion  may  be 
transmitted  from  the  rotating  shaft  to  the  flexible  casing  and  other  connecting  structures 
in  rotating  mechanical  equipment  [10-15],  These  simple  models  are  only  adequate  for 
the  free  and  forced  vibration  analyses  of  the  rotor  dynamic  system  enclosed  in  a  rigid 
casing.  For  example,  a  vibrational  model  of  a  rotating  system  based  upon  the  existing 
bearing  models  can  only  predict  purely  in-plane  type  motion  on  the  flexible  casing  plate, 
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given  only  the  bending  motion  on  the  shaft.  However,  experimental  results  have  shown 
that  the  casing  plate  motion  is  primarily  flexural  or  out-of-plane  type  [9,16,17].  This 
paradox  is  essentially  due  to  an  incomplete  understanding  of  the  bearing  as  vibratory 
motion  transmitter  in  rotating  mechanical  equipment 

The  main  focus  of  this  research  is  to  clarify  this  issue  quantitatively  and 
qualitatively  by  developing  a  new  mathematical  model  for  the  precision  rolling  element 
bearings,  and  extend  the  proposed  bearing  formulation  to  examine  vibration 
transmissibility  in  rotating  mechanical  equipment  through  several  example  cases  of 
bearing  systems  and  geared  drives.  The  superiority  of  the  proposed  model  compared  to 
simple  models  is  also  demonstrated  in  these  example  cases.  A  typical  shaft-bearing- 
casing-mount  system  is  shown  in  Figure  1.1.  The  rigid  or  flexible  shaft  may  be 
subjected  to  forces  and/or  torques  and  supported  by  a  bearing  on  a  flexibly  or  rigidly 
mounted  casing.  Here,  the  vibration  transmission  is  from  the  shaft  to  the  casing  and 
mount  through  the  bearing  system.  Figure  1.2  shows  a  typical  rolling  element  bearing 
subjected  to  forces  and  moments  due  to  the  rolling  element  deformation.  The  bearing  is 
free  to  rotate  about  the  axis  peipendicular  to  the  bearing  plane,  and  hence  does  not 
transmit  any  dynamic  moment  about  this  axis.  However,  dynamic  moments  about  the 
other  two  orthogonal  axis  exist  which  have  not  been  considered  in  simple  bearing 
models.  Finally,  a  generic  geared  rotor  system  consists  of  a  motor,  spur  gear  pair, 
flexible  shafts,  load,  flexible  couplings,  rolling  element  bearings,  flexible  casing  and 
mounts  as  shown  in  Figure  1.3a  is  also  considered.  The  system  is  excited  by  the  gear 
kinematic  transmission  error  at  gear  mesh  frequency.  A  single-stage  rotor  system  with 
rotating  mass  unbalance  excitation  as  shown  in  Figure  1.3b  is  treated  as  a  special  case  of 
Figure  1.3a.  Further  description  of  each  system  will  be  presented  in  later  chapters. 


Figure  1 . 1  Schematic  of  a  typical  shaft-bearing-plate-mount  system 
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Figure  1.2  Schematic  of  typical  rolling  element  bearings,  (a)  Bearing  forces  and 
moments,  (b)  Deep  groove  ball  bearing,  (c)  Angular  contact  ball  bearing, 
(d)  Taper  roller  bearing,  (e)  Cylindrical  roller  bearing. 
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Rigid  or  Flexible  Casing 


(a) 


(b) 


h 


Figure  1.3  Generic  rotating  mechanical  systems,  (a)  A  geared  rotor  system,  (b)  A 
single-stage  rotor  system 
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1.2  LITERATURE  REVIEW 

Simple  bearing  models  either  assume  ideal  boundary  conditions  for  the  shaft  or 
time-invariant  translational  springs  in  the  axial  and  radial  directions  [10-15].  The  ideal 
boundary  conditions  for  the  shaft  are  typically  simply-supported  for  short  bearings, 
clamped  for  long  bearings  or  free  (for  torsional  motion  only)  [10-12].  Formulas  for  the 
stiffness  coefficients  derived  from  the  force-displacement  relation  commonly  used  by 
bearing  designers  [18,19]  are  given  by  Harris  [18],  Gargiulo  [14]  and  White  [20].  In 
1982  Rajab  [21,22]  realized  the  limitations  of  the  simple  models  and  philosophically 
proposed  two  additional  stiffness  coefficients  which  couple  the  radial  and  rotational 
bearing  degrees  of  freedom,  given  radial  and  moment  about  the  axis  transverse  to  the 
radial  line  of  action.  In  1988,  Young  [23]  extended  Rajab’s  [21]  analyses  to 
include  mean  axial  force  while  retaining  other  features  of  Rajab’s  model.  This 
resulted  in  bearing  stiffness  matrix  of  dimension  3. 

Simple  bearing  models  are  widely  used  in  vibration  models  of  rotor  dynamic 
systems  which  typically  exclude  casing  and  mount  dynamics,  to  calculate  critical  speeds, 
dynamic  stability,  and  responses  due  to  shaft  excitations  such  as  mass  unbalance  and 
gear  transmission  error  [10-15].  In  most  cases,  the  vibration  transmission  through 
bearings  is  not  the  primary  issue,  and  thus  the  bearing  models  tend  to  be 
simplified.  None  of  the  current  models  studied  [20-22,24]  can  fully  explain 
vibration  transmission  through  bearings  in  systems  similar  to  Figure  1.1. 
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Publications  on  the  modal  analyses  of  geared  drives  and  single  or  multiple  stage 
rotor  systems  indicate  that  very  little  has  been  reported  directly  on  the  force 
transmissibility  through  bearings,  and  the  dynamic  effects  of  bearing,  casing  and 
mounts  on  the  internal  rotating  system  [10-13,15,25-29].  A  comprehensive  review  of 
the  available  vibration  studies  of  casing  and  mounts  has  been  given  by  Lim  and  Singh  in 
1989  [6].  In  these  studies,  the  dynamic  interaction  between  the  casing-mount  system 
and  the  internal  rotating  system  is  often  not  incorporated,  and  in  few  instances  when  this 
interaction  is  modeled,  only  purely  radial  and/or  axial  force  on  the  bearings  are  included. 
Such  models  still  do  not  explain  how  the  vibration  is  transmitted  from  the 
shaft  to  the  casing.  A  more  comprehensive  review  of  the  relevant  studies  will  be  given 
in  each  chapter. 

1.3  SCOPE  AND  OBJECTIVES 

A  new  mathematical  model  for  the  precision  rolling  element  bearing  in  Figure  1.2 
is  developed  and  incorporated  in  linear  time-invariant  discrete  and  broad  band  vibration 
models  of  Figures  1.1  and  1.3.  This  study  proposes  a  comprehensive  bearing  stiffness 
matrix  of  dimension  6  which  explains  the  vibratory  motion  transmission  through  the 
bearings  and  allows  for  the  study  of  overall  geared  rotor  system  dynamics.  The  lumped 
parameter  and  dynamic  finite  element  techniques  are  used  to  develop  the  discrete 
vibration  models  while  statistical  energy  analysis  method  is  used  for  the  broad  band 
vibration  models.  Experimental  validation  is  also  included  in  each  chapter,  where  the 
driving  point  and  cross  point  accelerance  and  mobility  levels  predicted  by  theory  are 
compared  to  experiments. 
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The  specific  objectives  of  this  research  are  grouped  as  follows:  bearing 
stiffness  formulation,  bearing  system  studies,  geared  rotor  system  studies,  and  statistical 
energy  analysis.  Each  chapter  is  self  sufficient  since  it  is  written  in  a  journal  paper  style. 
Accordingly,  a  detailed  problem  statement  is  also  included  in  each  chapter. 

a.  Bearing  stiffness  formulation:  Since  simple  bearing  models  used  in  rotor  dynamic 
analyses  are  inadequate  in  explaining  the  role  of  bearing  as  a  vibratory  motion 
transmitter,  this  study  resolves  this  issue  by  proposing  and  developing  a  new  rolling 
element  bearing  stiffness  matrix  which  is  suitable  for  the  analysis  of  the  vibration 
transmission  through  either  ball  or  roller  bearing.  A  numerical  scheme  is  also  developed 
to  compute  the  stiffness  coefficients  and  indicate  the  existence  of  solutions  to  the 
nonlinear  algebraic  bearing  equations  describing  the  bearing  load-displacement 
relationships.  The  proposed  bearing  model  is  partially  verified  by  comparing  with 
published  analytical  and  experimental  results.  In  addition,  the  character  of  the  bearing 
stiffness  matrix  and  its  sensitivity  to  various  bearing  parameters  will  be  discussed. 
(Chapter  II) 

b.  Bearing  system  studies:  The  specific  objectives  of  this  chapter  are  to  incorporate  the 
proposed  bearing  matrix  developed  in  Chapter  II  in  linear  discrete  vibration  models  of 
the  bearing  system  as  shown  in  Figure  1.1  using  lumped  parameter  and  dynamic  finite 
element  methods  to  compute  the  eigensolution  and  forced  harmonic  response,  and  to 
evaluate  the  dynamic  stability.  The  vibration  transmission  through  bearing  is  also 
predicted  for  several  example  cases  considered  previously  [14,20,24]  and  an 
experimental  setup  [17].  The  advantages  of  the  proposed  formulation  compared  to  the 
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simple  models  is  demonstrated  by  comparing  their  predicted  transfer  functions.  The 
theory  is  also  validated  by  comparing  analytical  predictions  with  experimental  data  on  a 
shaft-bearing-plate  system.  (Chapter  III) 

c.  Geared  rotor  system  studies:  Current  geared  rotor  system  vibration  models  exclude 
the  effects  of  casing  and  mounts,  and  do  not  address  the  overall  system  behavior.  The 
objectives  of  this  chapter  are  to  incorporate  the  proposed  bearing  matrix  in  the  discrete 
vibration  model  of  the  generic  geared  rotor  system  of  Figure  1.3  and  conduct  overall 
system  studies  by  calculating  eigensolutions  and  forced  harmonic  responses  with 
emphasis  on  the  prediction  of  vibration  transmission  through  rolling  element  bearings. 
The  effects  of  casing  and  mount  dynamics  on  the  internal  rotating  system  is  also 
evaluated.  Example  cases  which  include  a  single-stage  rotor  system  with  flexible  shafts 
supported  by  two  identical  rolling  element  bearings  on  rigid  casing  and  flexible 
mounts,and  a  spur  gear  pair  with  motor  and  load  inertias  attached  to  two  flexible  shafts 
supported  by  four  rolling  element  bearings  on  rigid  or  flexible  casing  and  compliant  or 
massive  mounts  will  be  studied  analytically  and/or  experimentally.  Also,  the  advantages 
of  the  proposed  formulation  as  compared  to  simple  models  of  geared  drives  will  be 
demonstrated.  (Chapter  IV) 

d.  Statistical  energy  analysis:  At  very  high  frequencies,  the  narrow  band  approach  using 
the  lumped  parameter  or  dynamic  finite  element  model  may  not  be  adequate  due  to  the 
high  structural  modal  density.  To  overcome  this  problem,  statistical  energy  analysis 
method  is  used  to  predict  the  vibratory  energy  transmission  in  and  noise  radiation  from  a 
geared  rotor  system  as  illustrated  in  Figure  1.4a.  The  proposed  bearing  matrix  is  again 
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Figure  1.4  Vibration  energy  transmission  in  a  geared  rotor  system,  (a)  Structure- borne 
noise  paths,  (b)  Statistical  energy  analysis  model. 
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incorporated  in  the  vibratory  energy  model.  In  this  method,  only  the  mean-square 
spatially  averaged  response  over  one  third  octave  frequency  bandwidths  are  predicted 
which  involves  solution  to  a  set  of  algebraic  equations  obtained  through  the  vibratory 
energy  balance  of  each  subsystem  shown  in  Figure  1.4b.  Several  example  cases 
including  a  plate-cantilevered  rectangular  beam,  circular  shaft-bearing-plate  system  and  a 
geared  rotor  system  are  chosen  to  demonstrate  the  salient  features  of  this  technique. 
(Chapter  V) 


CHAPTER  2 


BEARING  STIFFNESS  FORMULATION 

2.1  INTRODUCTION 

Current  rotor  dynamic  models  describe  precision  rolling  element  bearings  either  as 
ideal  boundary  conditions  for  the  shafts  [10-12],  or  as  purely  translational  stiffness 
elements  [13-15].  Such  simple  bearing  models  may  be  adequate  for  the  free  and  forced 
vibration  analyses  of  the  rotor  dynamic  system  enclosed  in  a  rigid  casing.  But  these 
mathematical  models  cannot  explain  how  the  vibratory  motion  may  be  transmitted  from 
the  rotating  shaft  to  the  flexible  or  rigid  casing  and  other  connecting  structures.  For 
example,  a  vibration  model  of  a  system  similar  to  Figure  2.1,  based  upon  the  existing 
bearing  models,  can  only  predict  purely  in-plane  type  motion  on  the  flexible  casing  plate 
given  only  the  bending  motion  on  the  shaft.  However,  experimental  results  have  shown 
that  the  casing  plate  motion  is  primarily  flexural  or  out-of-plane  type  [9,16,17].  This 
paradox  is  essentially  due  to  an  incomplete  understanding  of  the  bearing  as  vibratory 
motion  transmitter  in  rotating  mechanical  equipment  including  geared  drives  where 
structure-borne  noise  paths  through  bearings  are  often  dominant. 

This  chapter  clarifies  this  issue  qualitatively  and  quantitatively  by  developing  a 
new  mathematical  model  for  precision  rolling  element  bearings.  A  schematic  of  a 
generic  system  with  a  flexible  shaft  rotating  at  speed  Qz  and  subjected  to  mean  load 
vector  { f }  Sm=  ( Fwsm’Twsm ) ,  w  =  x,  y,  z,  flexible  casing  and  mount  is  shown  in  Figure 
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Figure  2.1  Schematic  representation  of  the  vibration  transmission  problem.  Here  the 
flexible  shaft  is  subjected  to  mean  forces  Fwsm  and  torques  Twsm  where  w 
=  x,  y  or  z,  is  the  direction  and  subscript  m  and  s  implies  mean  and  shaft 
respectively.  Also,  6  is  the  angular  displacement  and  u  is  the  translational 
displacement. 
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2. 1;  the  shaft  is  supported  on  one  of  the  following  bearings:  deep  groove  ball  bearing, 
angular  contact  ball  bearing,  thrust  ball  bearing,  straight  roller  bearing  or  taper  roller 
bearing.  A  new  bearing  stiffness  matrix  [K]bm  will  be  proposed  which  is  expected  to 
demonstrate  a  coupling  between  the  shaft  bending  motion  and  the  flexural  motion  of  the 
casing  plate.  It  will  be  shown  that  the  translational  bearing  stiffness  coefficients 
currently  used  in  rotor  dynamic  models  are  a  small  subset  of  the  proposed  [K]bm. 
Several  example  cases  are  employed  to  validate  our  theory.  Our  bearing  model  can  be 
easily  incorporated  in  analytical  or  numerical  models  typically  used  for  the  dynamic 
analyses  -  this  will  be  the  basis  of  Chapters  in  and  IV  of  this  report. 

2.2  LITERATURE  REVIEW 

The  ideal  boundary  conditions  for  the  shaft  have  typically  been  assumed  to  be 
simply-supported  for  short  bearings,  clamped  for  long  bearings  or  free  (in  the  torsional 
mode  only)  [10-12].  In  other  cases,  researchers  describe  the  bearing  as  time-invariant 
translational  springs  with  stiffness  coefficients  kbrr  and/or  kbzz  in  the  radial  and  axial 
directions,  respectively  [13-15,20].  Formulas  for  such  nonlinear  stiffness  coefficients 
are  given  by  Harris  [18]  and  Gargiulo  [14];  these  are  derived  from  the  radial  or  axial 
mean  force-displacement  equation  commonly  used  by  the  precision  rolling  element 
bearing  designers  [18,19].  Their  derivations  neglect  the  effects  of  radial  clearance  and 
mean  bearing  force  vector  [f)bm  on  the  load  distribution  and  hence  are  applicable  only 
for  constant  load  angle  y/(  of  180  degrees.  White  refined  these  formulations  by  using  a 
finite  difference  approximation  for  the  computation  of  stiffness  coefficients  for  radial 
ball  and  roller  bearings,  and  by  including  the  effects  of  radial  clearance  and  force  on  the 
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load  angle  \j//  [20].  Even  with  these  refinements  the  mathematical  model  is  still 
incapable  of  predicting  the  total  vibration  transmission  across  bearings. 

In  1982  Rajab  [21],  realized  the  limitations  of  the  current  simple  theory  and 
philosophically  proposed  two  additional  stiffness  terms  kbrQ  and  k^g  which  couple  the 
relative  radial  and  rotational  bearing  displacements  between  the  inner  and  outer  rings, 
given  the  mean  radial  load  and  moment  about  the  axis  transverse  to  the  radial  line  of 
action.  In  1988  Young  [23]  extended  Rajab's  [21]  analyses  to  include  the  mean 
axial  force  Fzbm,  and  then  used  a  discrete  summation  over  all  of  the  loaded 
rolling  elements  to  obtain  bearing  forces  and  moment  instead  of  the  integral  form 
while  still  retaining  other  features  of  Rajab's  model.  This  resulted  in  a  3x3 
bearing  stiffness  matrix.  Some  of  the  salient  features  of  Rajab's  [21]  and 
Young's  [23]  models  are  summarized  in  Reference  [22]. 

Experimental  determination  of  the  bearing  stiffness  coefficients  has  been  strictly 
limited  to  the  translational  coefficients  kbrr  and  kbzz.  A  method  for  the  measurement  of 
in  situ  bearing  stiffness  under  oscillating  loading  conditions  has  been  given  by  Walford 
and  Stone  [30].  Recently,  Kraus  et  al.  [24]  designed  an  in  situ  measurement  test  stand 
to  determine  the  translational  bearing  stiffness  from  measured  vibration  spectra,  in 
conjunction  with  the  single  degree  of  freedom  system  theory.  They  determined  the 
effect  of  preload,  bearing  release  and  rotational  speed  £2  z  on  kbrT  and  kbzz.  Their  results 
show  that  kbrr  and  kbzz  are  essentially  linear  and  the  effect  of  £2  z  is  negligible  when  a 
high  preload  is  applied  on  the  bearing. 
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2.3  ASSUMPTIONS  AND  OBJECTIVES 

Due  to  the  following  key  differences,  a  separate  formulation  of  [K]jjm  for  both  ball 
and  roller  type  rolling  element  bearings  is  required:  (i)  ball  bearings  have  elliptical 
contacts  and  roller  types  have  rectangular  contacts  between  the  inner  race,  rolling 
elements  and  outer  race  when  loaded,  and  (ii)  the  loaded  contact  angles  otj  of  the  ball 
types  may  change  but  otj  in  the  roller  type  remains  relatively  constant  [31].  Each  bearing 
is  characterized  by  its  kinematic  and  design  parameters  such  as  unloaded  contact  angle 
a0,  radial  clearance  rL,  effective  stiffness  coefficient  Kn  for  inner  ring-single  rolling 

element-outer  ring  contacts,  angular  misalignment,  preloads,  radius  of  inner  raceway 
groove  curvature  center  for  ball  type  and  bearing  pitch  radius  for  roller  type  [18,19,31]. 
It  is  expected  that  [K]bm  is  given  in  terms  of  these  parameters. 

The  mean  bearing  displacements  {q}bm  as  shown  in  Figure  2.2  are  given  by  the 
relative  rigid  body  motions  between  the  inner  and  outer  rings.  The  total  bearing 
displacement  vector  is  given  as  (q}b=(q}bm+{cl(t)}ba  where  [q(t))ba  »s  the  fluctuation 
about  the  mean  point  [q)bm  during  the  steady  state  rotation.  Accordingly  one  must 
consider  time  varying  bearing  stiffness  coefficients.  However  in  our  analysis,  such  time 
varying  bearing  stiffness  coefficients  are  neglected  by  assuming  very  small  vibratory 
motions  i.e.  [q)ba  «  (q)bm>  and  high  bearing  preloads.  Consequently,  only  the  mean 
bearing  loads  and  displacements  are  included  in  the  derivation  of  [KJbm.  The  basic 
load-deflection  relation  for  each  elastic  rolling  element  is  defined  by  the  Hertzian  contact 
stress  theory  [18,19,32],  and  the  load  experienced  by  each  rolling  element  is  described 
by  its  relative  location  in  the  bearing  raceway.  Further  it  is  assumed  that  the  angular 
position  of  each  rolling  element  relative  to  one  another  is  always  maintained  due  to  the 
rigid  cages  and  pin  retainers.  Secondary  effects  such  as  centrifugal  forces  and 


*  Outer  Ring  of  raceway  diamater 

*  Inner  Ring  of  raceway  diameter  d  bi 
Rolling  Element  of  diameter  D 


Figure  2.2  Rolling  element  bearing  kinematics  and  coordinate  system.  Here  the 
following  nomenclature  is  used:  dbo  is  the  outer  raceway  diameter,  dbm  is 
the  bearing  pitch  diameter,  is  the  inner  raceway  diameter,  y  is  the 
angular  position  of  rolling  element,  5wm  is  the  mean  translational 
displacement,  Ppm  is  the  mean  angular  displacement,  FWbm  is  the  mean 
bearing  force,  and  Mpbm  is  the  mean  bearing  moment  where  w  =  x,  y,  z, 
and  p=x,  y,  are  the  directions. 
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gyroscopic  moments  on  the  bearing  are  ignored  as  these  effects  are  evident  only  at 
extremely  high  rotational  speeds.  Tribological  issues  [32,33]  are  beyond  the  scope  of 
this  study  and  hence  our  analysis  assumes  bearings  to  be  unlubricated. 

The  specific  objectives  of  this  chapter  are  to:  (i)  propose  and  develop  a  new  rolling 
element  bearing  stiffness  matrix  [K]bm  which  is  suitable  for  the  analysis  of  the  vibration 
transmission  through  either  ball  or  roller  bearing,  (ii)  develop  a  numerical  scheme  to 
compute  [KJbn,  and  discuss  the  existence  of  solutions  to  the  nonlinear  algebraic  bearing 
equations  describing  load-displacement  relationships,  (iii)  verify  our  proposed  model  by 
comparing  its  predictions  with  published  analytical  and  experimental  results  [14,20,24] 
for  the  translational  stiffness  coefficients  kbxx,  kbyy  and  kbzz,  (iv)  relate  [K]bm  to 
various  kinematic  and  design  parameters,  and  perform  parametric  studies  to  investigate 
the  effect  of  unloaded  contact  angle  cXq  and  preloads,  and  (v)  characterize  the  nature  of 
[K]bm  and  recommend  its  usage.  Finally  it  should  be  noted  that  dimensionless 
parameters  will  not  be  used  here  as  the  metric  units  are  invariably  employed  to  specify 
bearings  [32], 

2,4  BEARING  LOAD-DISPLACEMENT  RELATIONS 

In  this  section,  the  relationships  between  the  bearing  forces  {Fxbm,  Fybm,  Fzbm) 
and  moments  (Mxbm,  Mybm}  transmitted  through  the  rolling  element  bearing,  and  the 
bearing  displacements  (q}bm  as  given  in  Figure  2.2  will  be  derived  for  both  ball  and 
roller  bearings.  The  mean  applied  loads  [fJsm  at  the  shaft  as  given  in  Figure  2.1  and 
bearing  preloads  generate  the  mean  bearing  displacements  {q)bm  and  loads  {f)bm- 
These  displacements  {q)bm  are  used  to  derive  the  resultant  elastic  deformation  Styj)  of 
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the  j-th  rolling  element  located  at  angle  Vj  from  the  x-axis.  From  the  ball  bearing 
kinematics  shown  in  Figure  2.3,  8g(\|/j)  is 


A(Vj)-A0  ,  8  >0 

5b(VJ)=|o  ,  5Bj£0 

A(Vj)  = 

(5*)zj  =  A0sin  a0+(8)zj;  (5*)jj  =  A0cos  a0  +  (8)rj 


(2.1a) 


(2.1b) 


(2.1c) 


where  A0  and  A  are  the  unloaded  and  loaded  relative  distances  between  the  inner  a;  and 
outer  &Q  raceway  groove  curvature  centers.  Similarly  for  the  roller  bearing  kinematics 
shown  in  Figure  2.4  for  aj=a0,  8R(Vj)  is 


(5)^  cos  a.  +  (5)  z.  sin  a. 
0 


6Rj>° 

8Rj£0 


(2.2) 


Note  that  in  equations  (2.1)  and  (2.2)  Sgj  £  0  or  5rj  :£  0  implies  that  the  j-th  rolling 
element  is  stress  free.  In  both  equations  (2.1)  and  (2.2),  the  effective  j-th  rolling 
element  displacements  in  the  axial  (6)zj  and  radial  (5)rj  directions  are  given  in  Figure  2.5 
in  terms  of  the  bearing  displacements  (qjbm- 

(5)  Zj  =  5zm  +  rj  {  Sin  j>  "  Pym  j) ]  <2‘3a> 

(5)rj=8xm  cos\|/j+  8ym  siny.-rL 


(2.3b) 
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Figure  2.3  Elastic  deformation  of  rolling  element  for  non-constant  contact  angle  ctj 
given  by  the  change  in  the  distance  between  the  inner  a;  and  outer  aQ 
raceway  groove  radius  curvature  centers  due  to  the  mean  bearing  loads  or 
displacements. 
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Figure  2.4  Elastic  deformation  of  rolling  element  for  constant  contact  angle  <Xj  =  a0 
given  by  the  change  in  the  relative  position  of  the  inner  and  outer  raceways 
due  to  the  mean  bearing  loads  or  displacements. 
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where  rj  is  the  radial  distance  of  the  inner  raceway  groove  curvature  center  for  the  ball 
type  or  is  the  pitch  bearing  radius  for  roller  type.  Equations  (2.1)-(2.3)  in  conjunction 
with  the  Hertzian  contact  stress  principle  [18,19,32]  stated  as  follows  yield  the  load- 
deflection  relationships  for  a  single  rolling  element . 

Qj  =  K„S”  (2.4) 

where  Qj  is  the  resultant  normal  load  on  the  rolling  element,  and  Kn  is  the  effective 
stiffness  constant  for  the  inner  race-rolling  element-outer  race  contacts  and  it  is  a 
function  of  the  bearing  geometry  and  material  properties  [18,19,31].  Note  that  the 
exponent  n  is  equal  to  3/2  for  ball  type  with  elliptical  contacts  and  10/9  for  roller  type 
with  rectangular  contacts.  Previously,  we  have  mentioned  that  the  loaded  contact  angle 
cxj  for  the  roller  bearing  remains  unchanged  from  the  unloaded  position  aG,  but  on  the 
other  hand  dj  may  alter  in  the  ball  bearing  case.  The  sign  convention  is  such  that  otj  is 
positive  when  measured  from  the  bearing  x-y  plane  towards  the  axial  z-axis  as  shown  in 
Figures  2.3  and  2.4,  and  negative  otherwise.  For  the  ball  bearing  of  Figure  2.3,  the 
loaded  contact  angle  Otj  is 


tan  (a.)  = 


A0  sin  a0  +  (8)zj 
A  o  cos  a0  +  (8)^ 


(2.5) 


where  (5)zj  and  (8)rj  are  given  by  equations  (2.3a)  and  (2.3b).  It  is  appropriate  here  to 
note  that  Rajab  [21,22]  and  Young  [22,23]  in  their  derivation  of  the  bearing  stiffness 
model  used  an  expression  similar  to  equation  (2.2)  but  with  8xm=Pym=8zm=rL=0  in 
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Rajab's  analysis  and  5xrn=Pyrn=rL=0  in  Young's  analysis  for  both  ball  and  roller 
bearings.  Since  always  otj  is  given  by  equation  (2.5)  irrespective  of  the  formulation  and 
since  equation  (2.2)  is  valid  only  if  ctj=aQ,  their  ball  and  roller  bearings  analyses  are 

in  error.  Expressions  similar  to  equation  (2.1)-(2.5)  with  minor  differences 
have  also  been  used  by  Eschmann  et  al.  [31],  Jones  [34]  and  Davis  [35],  but  their 
intentions  were  to  calculate  static  bearing  forces  rather  than  to  derive  the  bearing  stiffness 
models  for  vibration  transmission  analysis. 


2.5  DEVELOPMENT  OF  BEARING  STIFFNESS  MATRIX  [K]bm 

Our  proposed  bearing  stiffness  matrix  [K]bm  is  a  global  representation  of  the 
bearing  kinematic  and  elastic  characteristics  as  it  combines  the  effects  of  z  number  of 
loaded  rolling  element  stiffnesses  in  parallel  given  by  5j  >  0.  First,  we  need  to  relate  the 
resultant  bearing  mean  load  vector  { f }bm  t0  t*ie  bearing  displacement  vector  [q}bm.  This 
can  be  achieved  through  vectorial  sums  Qj  (8wm,Ppm;  w  =  x,  y,  z  and  p  =  x,  y)  in 
equation  (2.4)  for  all  of  the  loaded  rolling  elements  which  lead  to  the  following  bearing 
moments  j  Mwbfn }  and  forces  { Fwbm }  as  follows 


M  h 
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sin  Vj 

M  . 

ybm 

•  =  y,  r .  Q  .  sin  a .  ■ 
“  J  J  ) 

J 

-  cos  \|T 

M  . 

zbm 

0  . 

xbm 

z 

cos  a .  cos  \ir . 

J  TJ 

ybm 

cos  a.  sin  w. 

J  YJ 

zbm 

sin  a. 

J 

(2.6a) 


(2.6b) 
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Replacing  Qj  and  otj  in  equation  (2.6)  in  terms  of  {8wm,Ppm}  yields  the  following 
explicit  relationships  between  {f}^  and  (q)bm  f°r  ball  bearings 


u 


[AQ  sin  a0+(8)zj]2+[A, 


_  n 

c°sa0+(8)i.]2-A0} 


[A0sina0+(8)zjr 


+  [A0  cosa0+(5)|;.] 


r.  {A  0  sin  a0+  (8)zj) 


f 


-  cos  y  • 


0  J 


(2.7a) 


j- [A0  cos  a0  +  (8)r.]  cos  \|r  ' 

.  [A0  cos  a0  +  (8)r.]  sin  y.  . 
[A0sin  a0  +  (S)zj] 


(2.7b) 


and  similarly  for  roller  bearings 
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where  (5)rj  and  (8)zj  are  functions  of  {5wm,ppm}  as  defined  by  equation  (2.3). 
Approximate  integral  forms  of  equations  (2.7a,b)  and  (2.8a,b)  are  often  used  instead  of 
the  summation  forms  to  eliminate  explicit  dependence  on  yj,  especially  in  the  case  of 
only  one  or  two  degrees  of  freedom  bearings  [18,32].  For  instance  Rajab  [21,22]  chose 
the  integral  form  representation  but  made  a  mathematical  error  in  constructing  the 
integrand. 

Now  we  define  a  symmetric  bearing  stiffness  matrix  [K]bm  of  dimension  6  from 
equations  (2.7a,b)  and  (2.8a,b)  and  by  assuming  that  [q)ba  «  {q)bm 


w,  i  =  x,  y,  z  (2.9) 


Here  each  stiffness  coefficient  must  be  evaluated  at  the  mean  point  [qJbm*  Explicit 
expressions  for  the  ball  bearing  stiffness  are  as  follows;  note  that  [K]bm  is  symmetric 
i.e.  kbiw=kbwi- 
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where  (8  )zj,  (8  )rj  and  Aj  are  defined  by  equation  (2.1).  And  the  roller  bearing 
stiffness  coefficients  kbjw  =  kbwj  are  given  explicitly  as 
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where  8Rj  is  defined  in  equation  (2.2).  It  should  be  noted  that  all  stiffness  terms 
associated  with  the  torsional  degree  of  freedom  are  zero  due  to  the  fact  that  an  ideal 
bearing  allows  free  rotation  about  the  z-direction.  Also,  the  translational  stiffness 
coefficients  kbii.  i=x,y,z  for  6ym=8zm=pxm=pym=0,  5xm=5zm=Pxm=pym=0  or 
8xm=8ym=Pxm=Pym=^  are  equivalent  to  the  bearing  stiffness  coefficients  commonly 
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used  by  investigators  [14,15,20].  The  nature  of  these  and  other  features  of  [K]bm  wid 
be  discussed  later  in  Section  2.9. 

2.6  NUMERICAL  ESTIMATION  OF  [K]bm 

The  coefficients  kbiw  can  be  computed  by  one  of  the  following  two  methods:  I. 
directly  compute  kbiw  given  mean  bearing  displacement  vector  {q}bm  employing 
equations  (2.10a-p)  and  (2.11a-p),  or  II.  numerically  solve  the  nonlinear  algebraic 
equations  described  by  equations  (2.7a, b)  and  (2.8a,b)  to  obtain  [q}bm  from  {f}bm’ 
and  then  evaluate  kbiw  per  method  I.  Note  that  {f)bm  may  be  functions  of  the  mean 
shaft  loads,  bearing  preloads,  and  shaft  and  casing  compliances  depending  on  the 
configuration  and  flexibility  of  the  rotating  mechanical  system.  If  the  bearing  system  is 
statically  determinate,  then  {f}bm  may  be  computed  explicitly  in  terms  of  [f}sm  and 
preloads  using  the  force  and  moment  equilibrium  equations.  Conversely  for  an 
indeterminate  system,  appropriate  field  equations  for  the  shaft  and  casing  plate  are 
needed  in  addition  to  the  equilibrium  equations  to  obtain  {fJbm  which  must  also  include 
shaft  and  casing  compliances.  Calculations  of  [f}bm  (QJbm  *n  case  are 
simultaneous,  which  may  be  extensive  especially  when  the  system  is  very  flexible,  and 
may  even  require  discretization  using  finite  element  or  lumped  mass  technique. 
However,  in  many  real  machines  the  in-plane  stiffness  of  the  casing  plate  which 
supports  most  of  the  mean  bearing  load  is  much  higher  than  the  bending  stiffness  of  the 
shaft.  Hence  the  casing  in-plane  stiffness  term  may  be  neglected  without  contributing 
any  large  error  to  {f}bm  [18,31].  And  only  the  Euler’s  beam  equation  for  a  statically 
indeterminate  shaft  is  used  along  with  the  nonlinear  bearing  load-displacement  equations 
(2.7a, b)  and  (2.8a,b). 
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Method  I  is  computationally  direct  and  needs  no  discussion.  But  method  II  deals 
with  as  many  as  10  N  nonlinear  algebraic  equations  for  N  bearings  if  the  casing 
flexibility  is  neglected.  One  must  choose  an  appropriate  numerical  method  as  the 
nonlinear  algebraic  equations  must  be  solved  iteratively  [36,37].  In  addition,  the 
available  numerical  methods  need  a  prior  knowledge  of  the  approximate  location  of  the 
solution  vector  being  sought  and  hence  one  must  be  careful  in  interpreting  the  numerical 
results.  In  this  study,  we  adopted  the  Newton-Raphson  method  for  its  good 
convergence  characteristic  [36,37].  To  implement  this  method,  equation  (2.6)  for  each 
bearing  is  rearranged  as 
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where  Hi,H2,...,H5  are  functions  defined  for  computational  reasons.  For  an 

indeterminate  system,  there  are  additional  functions  H6,H7,...,Hv  from  the  field 

equations.  Using  Taylor’s  series,  any  function  in  equations  (2.12a,b)  can  be 

expanded  about  the  solution  vector  X  =  [q]bm  for  a  statically  determinate  system  andX 
T  t  T 

=  ^  ^  bm’  ^  bm  ^  a  stat‘cally  indeterminate  system  as  follows  by  neglecting 


second  and  higher  order  terms. 
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V 

Hk  (X+SX)  -  Hk  (X)  +  £  8x .  ;  k=l,2,3,...V  (2.13) 

j  i 

The  solution  for  the  incremental  vector  8X  can  be  obtained  by  setting  Hk(X+8X)  =  0 
per  equations  (2.12)  and  (2.13)  which  yields  a  set  of  linear  algebraic  equations.  This 
vector  8X  is  added  to  the  previously  computed  vector  X  given  by  Hk(X)  =  0  for  the 
next  iteration  until  the  convergence  criterion,  say  that  8X  is  within  a  specified  tolerance, 
is  satisfied.  Our  proposed  numerical  scheme  can  be  summarized  as  follows:  (i)  guess 
bearing  displacement  vector  {q}bm  and/or  load  vector  {f)bm»  (“)  compute  8X  and  check 
against  a  specified  tolerance,  (iii)  add  8X  to  the  previous  solution  vector  X  and  repeat 
steps  (i)  and  (ii)  until  the  convergence  criterion  is  satisfied.  We  have  found  that  a  few 
initial  guess  trials  are  required  in  most  cases  to  obtain  reasonable  results. 

2.7  VALIDATION  OF  PROPOSED  MODEL 

In  order  to  validate  our  theory  we  compare  the  translational  stiffness  coefficients  of 
the  proposed  bearing  matrix  [K]bm  with  published  analytical  and  experimental  results 
[14,20,24].  First  we  apply  our  theory  to  predict  the  nonlinear  axial  kbzz  =  kbzz(8zm) 
and  radial  kbrr  =  kbrr(8rm)  stiffnesses  as  shown  in  Figure  2.6.  Our  predictions  are 
found  to  be  within  2%  of  Gargiulo's  [14]  formulas  which  are  commonly  used  for  both 
ball  and  roller  bearings. 

For  the  second  example  case,  we  consider  the  ball  bearings  used  by  Kraus  et 
al.[24]  for  an  in-situ  determination  of  the  bearing  stiffness.  Using  their  bearing  design 
parameters,  we  compute  radial  stiffness  coefficient  k^  as  a  function  of  the  axial  preload 
Fzbm.  Excellent  comparison  between  theory  and  experiment  is  seen  in  Figure  2.7. 


Radial  Stiffness  (MN/mm) 


Figure  2.6  Comparison  between  the  proposed  theory  and  Gargiulo's  formulas  [14]  for 
axial  kbzz  and  radial  stiffness  coefficients  of  ball  and  roller  bearings. 
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Axial  Preload  (kg) 


Figure  2.7  Comparison  between  the  proposed  theory  and  the  experimental  results  of 
Kraus  et  al.[24]  for  k^  as  a  function  of  the  mean  axial  preload. 
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Finally,  we  compare  our  results  for  the  nonlinear  radial  stiffness  kblT  with  those 
reported  earlier  by  White  [20]  for  both  ball  and  roller  bearings.  We  note  discrepancies 
in  Figure  2.8  between  our  theory  and  White's  results.  In  order  to  explain  these  we  now 
define  kbrr  using  the  finite  difference  approximation  which  was  also  used  by  White: 
k  brr  ~  AFbrm  /  ASm,  =  Fbrm  /  (8rm-rL).  Now  a  good  match  is  evident  in  Figure  2.8 
between  our  kbrr  values  and  the  data  given  by  White.  However,  the  correct  formulation 
is  obviously  given  by  our  proposed  theory  which  is  based  on  the  analytical  partial 
derivatives  k^  =  3Fbrm  /  3801,  as  the  displacement  8m  may  be  large. 

2.8  PARAMETRIC  STUDIES 

The  proposed  matrix  [K]bm  includes  a  coupling  between  the  casing  flexural 

motion  and  shaft  bending  motion  which  is  reflected  by  some  of  the  dominant  off- 
diagonal,  kbx0y,  kbyex,  kbz9x  and  kbz0y,  and  rotational  diagonal,  kb0x0x  and  kb0y0y, 

stiffness  coefficients;  these  are  labeled  as  'coupling  coefficients'  for  discussion 
purposes.  Such  stiffness  coefficients  are  investigated  further  by  varying  preloading 
conditions  and  unloaded  contact  angle  Oq  for  both  ball  (set  A)  and  roller  (set  B)  bearings 
whose  design  data  are  listed  in  Table  2.1. 

The  coupling  coefficients  given  a  constant  mean  radial  displacement  8rm  (radial 
preload),  as  shown  in  Figures  2.9  and  2.10  for  both  ball  and  roller  bearings 
respectively,  are  found  to  increase  as  a0  increases  and  reach  a  maximum  when  0^  is 
near  90°.  On  the  other  hand,  the  radial  translational  stiffness  coefficients  in  the  x  and  y 
directions  are  found  to  decrease  as  Oq  increases.  These  observations  imply  that  for  deep 
groove  ball  type  or  straight  roller  type  bearing  (Oq  =  0°)  the  radial  stiffness  coefficients 
kbrr  are  dominant,  but  for  angular  contact  ball  type  or  taper  roller  type  bearing  (Oq  >  0°) 


Radial  Stiffness  (kN/mm) 


200  y 


Proposed  Theory 

Estimated 

White 


Radial  Deflection  (mm) 


Figure  2.8  Comparison  between  the  proposed  theory  kj„p  estimated  k^,  and  White’s 
analytical  results  [20].  a.  Ball  bearing,  b.  Roller  bearing. 


Parameters 


Set  A  (ball  type) 


Set  B  (roller  type) 


Load-deflection  exponent  n 

3/2 

10/9 

Load-deflection  constant  Kn  (N/mn) 

8.5  E9 

3.0  E8 

Number  of  rolling  element  Z 

12 

14 

Radial  clearance  rL(mm) 

0.00005 

0.00175 

Pitch  radius  ^(mm) 

19.65 

21.25 

Ao  (mm)  * 

0.05 

- 

t  Unloaded  distance  between  inner  and  outer  raceway  groove  curvature  centers  (see 
Figure  2.3) 

tt  Equivalent  to  rj  for  roller  bearings  and  rj-Ao/2  for  ball  bearings  given  in  equation 
(2.3) 


Bearing  Stiffness  (MN/rad)  Bearing  Stiffness  (kNm/rad)  Bearing  Stiffness  (kN/mm) 
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Figure  2.9  Dominant  stiffness  coefficients  of  ball  bearing  set  A  for  0°  £  Oq  £  90°  and 
given  a  constant  mean  radial  bearing  displacement  S^, ■  0.025  mm. 


Bearing  Stiffness  (GN/rad)  Bearing  Stiffn< 


Figure  2.10  Dominant  stiffness  coefficients  of  roller  bearing  set  B  for  0°  £  a0  2  90° 
and  given  a  constant  mean  radial  bearing  displacement  Sxm  =  0.025  mm. 
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the  coupling  terms  are  more  significant.  Note  that  in  Figure  2.10,  all  the  stiffness 
coefficients  are  zero  at  ac  =  90°  for  the  roller  type.  This  is  due  to  the  fact  that  in  the 
thrust  roller  bearing,  radial  flanges  are  included  to  resist  the  roller  motion  in  this 
direction  which  is  not  modeled  here,  and  hence  these  stiffness  coefficients  must  vanish. 
In  addition,  thrust  roller  bearings  are  designed  to  carry  axial  loads  [18,31].  On  the  other 
hand,  ball  bearings  have  finite  stiffness  coefficients  at  Oq  =  90°  due  to  the  curvature  of 
the  raceway  which  provide  some  resistance  to  the  radial  preloads.  In  general,  the  trends 
in  both  ball  and  roller  bearing  stiffness  properties  are  similar  when  each  is  subjected  to 
mean  radial  displacement  or  preload. 

In  the  case  when  the  bearings  are  subjected  to  mean  axial  displacement  (axial 
preload),  as  shown  in  Figure  2.11  for  the  ball  type  and  Figure  2.12  for  the  roller  type, 
the  number  of  nonzero  stiffness  coefficients  are  less  than  those  seen  for  the  radial 
preload  only.  Again,  it  is  observed  that  both  ball  and  roller  bearings  display  similar 
trends.  Over  mid  to  high  a0  values,  the  coupling  coefficients  are  found  to  be 
significant.  The  translational  stiffness  coefficients  are  relatively  constant  except  for  the 
axial  stiffness  which  increases  as  (Xq  increases.  This  is  expected  due  to  the  inclination  of 
the  rolling  element  line  of  contact  from  the  x-y  plane  which  increases  elastic  support  in 
the  z-direction.  At  cLq  =  0°,  all  the  stiffness  coefficients  for  roller  bearings  are  zero  as 
there  is  no  constraint  in  the  axial  direction.  In  real  bearings  such  a  constraint  is  provided 
by  the  axial  flanges  [18,31],  however  this  bearing  is  not  designed  to  carry  any  axial 
preload. 

Results  for  the  misalignment  in  ball  and  roller  bearings  simulated  by  specifying  a 
mean  bearing  angular  displacement  Pym  are  shown  in  Figures  2.13  and  2.14 
respectively.  The  dominant  stiffness  coefficients  are  the  same  as  those  seen  for  the 


Bearing  Stiffness  (GN/rad)  Bearing  Stiffness  (GNm/rad)  Bearing  Stiffness  (MN/mm) 


Bearing  stiffness  (MN/rad)  Bearing  Stiffness  (MNm/rad)  Bearing  Stiffness  (MN/mm) 
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Unloaded  Contact  Angle  (deg) 


Figure  2. 1 3  Dominant  stiffness  coefficients  of  ball  bearing  set  A  for  0°  £  (Xq  £  90°  and 
given  a  constant  misalignment  Pym  =  0.015  rad. 
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radial  preload  case.  For  ball  bearings,  most  of  the  stiffness  coefficients  remain  constant 
for  0°  <  Oq  <  90°.  On  the  other  hand,  the  stiffness  coefficients  for  roller  bearing  have 
trends  similar  to  those  found  for  the  radial  preload  cases. 

From  the  detailed  parametric  studies,  it  is  concluded  that  the  nature  of  [K]bm  is 
dictated  by  the  bearing  type,  Oq  and  preloads.  Also,  the  coupling  coefficients  are  not 
negligible  in  most  cases  as  assumed  previously  by  many  investigators. 

2.9  CONCLUSIONS 

Results  of  Section  2.8,  which  show  similar  trends  for  some  of  the  cases,  imply 
that  there  may  be  a  systematic  approach  to  characterize  the  proposed  bearing  stiffness 
matrix  [K]bm.  From  the  kinematic  and  geometrical  considerations,  it  is  always  possible 
to  impose  any  bearing  displacement  vector  {q}bm  which  denotes  relative  rigid  body 
motions  between  the  inner  and  outer  rings  as  long  as  the  rolling  element  is  still  within 
the  elastic  deformation  regime.  On  the  other  hand,  ah  arbitrary  application  of  {f}bm  may 
not  produce  a  singular  displacement  response  from  the  bearing  due  to  its  kinematic  and 
geometrical  constraints.  Hence,  we  compute  [K]bm  and  {f}bm  by  systematically 
varying  {q}bm.  The  results  of  all  possible  forms  of  [K]bm  are  listed  in  Table  2.2  and 
2.3  for  ball  and  roller  bearings  respectively.  Also  included  here  are  the  current  bearing 
models  which  are  based  on  the  translational  spring  descriptions;  these  models  do  not 
show  any  coupling.  Note  that  the  exact  values  of  the  stiffness  coefficients  are  not  given 
as  these  depend  on  specific  parameters;  therefore  only  the  dominant  kyj  terms  are  listed 
for  all  possible  bearing  load  configurations  along  with  the  corresponding  {q}^  and  Oq. 
Also,  note  that  not  all  combinations  of  the  bearing  loads  are  possible  which  complicates 
bearing  stiffness  calculations  further,  especially  for  the  numerical  method  II.  Tables  2.2 


Table  2.2  Comparison  between  the  proposed  and  current  ball  bearing  stiffness 
coefficients.(p  =  x,  y;  i  =  x,  y  but  i  *  p) 


Mean 

Mean  bearing  displacement 

Dominant  stiffness  coefficients 

bearing^ 

loads 

«„"0° 

0°<ao<90° 

ao«90° 

current^ 

proposed^* 

Ppm 

^pm 

- 

- 

kPP 

^xx»^yy»^zz>^0x0x’^0y0y’^z0i 

Pzm 

^zm 

^zm 

— 

^pp»^ZZ 

^xx»^yy»^zz>^0x0x,k0y0y’ 

^X0y>^y0x 

Fzm 

- 

- 

^zm 

kzz 

^xx^yy^zz’^OxOx’^OyOy 

Mpm 

Ppm 

- 

- 

- 

kxx^yy^zz^OxOx^OyOy’kiz 

Pzm*^pm 

- 

- 

^zm’Ppm 

^zz 

^xx^yy’^zz'^OxOx’^OyOy’^zOp 

Pxm’Pym 

®xm>®yra 

— 

— 

kpp 

^xx>^yy>^zz>^0x0x’^0y0y’ 

^xy’^xOy’^zOx^zOy 

F  M 

1  pm>*vlpm 

^pm>Ppm 

““ 

kPP 

^xx*^yy>^zz’^0x0x,k0y0y’ 

kx0x’ky0y’kpz>kz0i 

^pm>Ppm 

— 

“ 

^xx>^yy>^zz>^0x0x,k0y0y* 

^xy^xz’^yz.^OxOy 

Ppm»Pzm* 

^im 

^pm>^zm> 

Pirn 

Spm>^zm> 

Pim 

Spm’^zm* 

Pim 

lr  Jr 
^ppj^zz 

^xx»^yy>^zz»^0x0x,k0y0y* 

^x0y»^y0x,kpz,^z0i 

Pzm>Mxm> 

^lym 

^zm>Pxm* 

Pym 

®zm»Pxm> 

Pym 

^zm*Pxm* 

Pym 

^PP’^ZZ 

^xx>^yy*^zz’^0x0x,k0yOy’ 

^xy>^0x6y,kz0x*^z0y 

combinations  of  {q}m 

kppt^ZZ 

all  non-zero  except  6Z  terms 

t  Ideal  boundary  condition  models  used  to  describe  the  bearing  are  not  tabulated. 

tt  Here  the  subscript  b  which  implies  bearing  has  been  omitted  for  brevity.  j 

ttt  All  terms  associated  with  0Z  are  zero  because  of  the  free  rotation  about  the  z  axis.  $ 

i 

;  C 
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Table  2.3  Comparison  between  the  proposed  and  current  roller  bearing  stiffness 
coeflficients.(p  =  x,  y;  i  =  x,  y  but  i  9*  p) 


Mean 

bearing^ 

loads 

Mean  bearing  displacement 

Dominant  stiffness  coefficients 

s 

°« 

0 

0 

0°<ao<90° 

ao«90° 

current^ 

proposed^ 

Fpm 

5pm 

- 

- 

kPP 

^bpp 

Fzm 

— 

5zm 

— 

kpp^zz 

^xx»^yy*^zz>^0x0x,k6y0y’ 

^x0y>^y0x 

F 

rzm 

- 

- 

5zm 

^zz 

^zz>^6x0x,k0y0y 

Fzm’Mpm 

- 

- 

5zm>Ppm 

^zz 

^zz>^0x0x,k0y0y,kzOp 

Fxm>Fym 

5xm»5ym 

- 

- 

kpp 

^xx»^yy»^xy 

Fpm>Fzm 

5pm>5zm 

Pim 

*** 

kpp»kzz 

^xx>^yy>^zz>^0x0x'^0y0y’ 

^X0y»^y0x,kpz,kz0i 

F  M 

1  zm*iV1xm» 

Mym 

“ 

5zm>Pxm> 

Pym 

Kz 

^zz»^0x0x,kOy0y,k0x0y’ 

^Z0x,kz0y 

(On, 

combinations  of  {q}m 

kpp4tzz 

all  non-zero  except  0Z  terms 

t  Ideal  boundary  condition  models  used  to  describe  the  bearing  are  not  tabulated. 

tt  Here  the  subscript  b  which  implies  bearing  has  been  omitted  for  brevity. 

ttt  All  terms  associated  with  0Z  are  zero  because  of  the  free  rotation  about  the  z  axis. 
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and  2.3  should  provide  some  insight  to  the  solution  of  the  nonlinear  algebraic  bearing 
load-deflection  equations  which  requires  a  prior  knowledge  of  the  type  of  solution 
being  sought  as  outlined  earlier.  In  most  practical  problems,  mean  bearing  loads  are 
typically  known.  This  knowledge  can  be  combined  with  Table  2.2  or  2.3  to  formulate 
the  nonlinear  load-deflection  equations  in  the  simplest  form  by  deleting  all  of  the  zero 
displacement  terms. 

Tables  2.2  and  2.3  show  that  the  coupling  coefficients  k|,x0y,  kby0x,  kbz0x,  kbz0y, 
kb0x0x  and  are  found  to  be  dominant  in  most  of  the  ball  bearing  cases,  and  only 

in  some  of  the  roller  bearing  cases.  This  is  essentially  due  to  the  curvature  of  the 
raceway  in  ball  bearing  which  invariably  causes  the  rolling  element  to  orient  itself  such 
that  0°  <  otj  <  90°  which  generates  ball  loads  in  the  z  direction  as  well.  However,  in  the 
roller  bearing  case  where  otj  =  a0,  the  same  phenomenon  does  not  occur  when  a0  =  0° 
or  90°,  and  the  coupling  coefficients  are  seen  only  when  ctj,  *  0°  or  90°.  In  fact  for  the 
0°  and  90°  unloaded  contact  angle  cases,  the  stiffness  coefficients  associated  with  x  and 
y  directions  and  those  associated  with  the  z,  0X,  and  0y  directions  do  not  exist 
simultaneously;  the  former  is  dominant  when  Oq  =  0  and  the  latter  prevails  when  oCq  = 
90°.  Another  case  of  interest  here  is  the  case  when  bearing  loads  are  complex  as  given 
by  the  last  row  in  Tables  2.2  and  2.3  where  all  of  the  bearing  stiffness  coefficients 
unrelated  to  the  rotational  degree  of  freedom  0Z  exist.  Solution  to  these  cases  may 
require  a  large  number  of  iterations. 

In  summary,  we  have  developed  a  comprehensive  bearing  stiffness  matrix  from 
the  basic  principles  which  includes  all  possible  rigid  body  degrees  of  freedom  of  a 
bearing  system.  This  matrix  has  been  validated  partially  using  several  analytical  and 
experimental  examples.  Further  validation  of  [KJbm  is  not  possible  as  coupling 
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coefficients  are  never  measured  [24,30].  Nonetheless,  our  theory  is  general  in  nature 
and  is  applicable  to  even  those  configurations  which  may  be  different  from  the  generic 
case  shown  in  Figure  2.1.  Further  research  is  required  to  incorporate  tribological  issues 
[32,33]  in  this  formulation.  However  the  proposed  stiffness  matrix  in  its  present  form, 
unlike  the  current  models,  is  clearly  capable  of  explaining  the  nature  of  vibration 
transmission  through  bearings  -  this  is  the  subject  of  Chapters  III  and  IV  of  this 
report,  which  will  also  include  further  comparisons  between  theory  and  experiment. 


CHAPTER  HI 


BEARING  SYSTEM  STUDIES 

3.1  INTRODUCTION 

Current  bearing  models  [10-15]  can  not  explain  how  the  vibratory  motion  may  be 
transmitted  from  the  rotating  shaft  to  the  casing  and  other  connecting  structures  in 
rotating  mechanical  equipment.  For  instance,  experimental  results  [9,16,17]  have 
shown  that  casing  plate  motion  for  a  system  similar  to  Figure  3.1  is  primarily  flexural  or 
out-of-plane  type  given  only  the  bending  motion  on  the  shaft  Using  existing  vibration 
models,  only  in-plane  type  motions  on  the  casing  plate  are  obtained.  Such  limitations 
associated  with  current  bearing  models  have  been  discussed  thoroughly  in  Chapter  II  of 
this  report  Also  in  Chapter  II,  a  new  mathematical  model  for  the  precision  rolling 
element  bearings  has  been  developed  in  order  to  clarify  this  issue  qualitatively  and 
quantitatively. 

This  study  extends  the  proposed  bearing  formulation  and  demonstrates  its 
superiority  over  the  existing  models  in  vibration  transmission  analyses.  A  schematic  of 
a  generic  system  with  a  flexible  shaft  rotating  at  constant  speed  Cl  z,  flexible  casing  and 

mount  is  shown  in  Figure  3.1.  The  shaft  is  supported  by  a  rolling  element  bearing 
which  is  modeled  by  a  stiffness  matrix  [K]j,m  of  dimension  6  as  proposed  in  Chapter  n. 
The  excitations  at  the  rotating  shaft  are  given  in  terms  of  an  alternating  load  vector 
(f(t))sa  *  {FjsafO.Tjsaft))7  *  {f(t))g-  {fjgn,;  j«x,y,z,  where  FjsaCt)  and  Tj^ft)  are  the 


52 


Figure  3.1  Schematic  representation  of  the  vibration  transmission  problem.  Here  the 
flexible  shaft  is  subjected  to  alternating  forces  Fj^t)  and  torques  TjsaW 
where  j  =  x,  y  or  z,  is  the  direction  and  subscript  a  implies  alternating. 
Also,  6  is  the  angular  displacement  and  u  is  the  translational  displacement 
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alternating  force  and  torque  respectively,  {f(t)}s  is  the  total  load  vector  of  dimension  6, 

{ f }  sm  represents  the  mean  load  vector,  and  superscript  T  implies  the  transpose.  In  the 
vibration  analysis,  (flsm  and  bearing  preloads  are  not  included  as  they  do  not  appear  in 
the  governing  equations  of  the  linear  vibration  model  but  are  used  for  computing  [K]t,m. 
The  effect  of  bearing  coupling  coefficients,  which  are  off-diagonal  and  rotational 
diagonal  terms  of  [K]f,m  as  described  in  Chapter  II,  on  the  eigensolution,  forced 
vibration,  and  vibration  transmission  through  bearings  is  evaluated.  Our  theory  will  be 
illustrated  and  validated  through  3  physical  system  example  cases;  experimental 
verification  is  also  included. 

3.2  LITERATURE  REVIEW 

The  existing  bearing  models  which  assume  either  ideal  boundary  conditions  [10- 
12]  for  the  shaft  or  translational  stiffness  elements  [13-15]  have  already  been  discussed 
in  Chapter  II.  Various  formulas  for  estimating  translational  stiffness  coefficients 
commonly  used  by  researchers  have  been  compared  with  our  proposed  [K]bm 
formulation.  These  simple  bearing  models  are  widely  used  in  vibration  models  of  the 
rotor  dynamic  systems,  which  typically  exclude  casing  and  mount  dynamics,  to  calculate 
critical  speeds,  responses  due  to  shaft  excitations  such  as  mass  unbalance  and  gear 
transmission  error,  and  dynamic  stability  [10-15].  In  most  of  these  cases,  the  vibration 
transmission  through  bearings  is  never  or  not  the  primary  issue,  and  thus  the 
bearing  models  tend  to  be  simplified.  None  of  the  current  models  studied  [20-22,24] 
can  fully  explain  vibration  transmission  through  bearings  in  systems  similar  to 
Figure  3.1.  In  1979  White  [20]  evaluated  the  rolling  element  bearing  vibration  transfer 
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characteristics  using  a  two  degrees  of  freedom  (DOF)  vibration  model  of  the  system 
shown  in  Figure  3.1.  His  formulation  is  based  on  only  the  radial  bearing  stiffness 
coefficient  kbrr .  He  concluded  that  an  increase  in  preload  increases  kbrr  and  system 
natural  frequencies.  He  also  found  that  the  effect  of  bearing  nonlinearity  is  negligible  at 
higher  preloads.  In  1987  Kraus  et  al.  [24]  proposed  a  single  degree  of  freedom  model 
for  a  similar  physical  system  (with  a  very  compliant  mount)  to  estimate  kbrr  from 
measured  vibration  transmission  spectra.  In  both  of  these  studies,  the  coupling 
coefficients  of  [K]bm  arc  not  included. 

In  1982  Rajab  [21]  philosophically  proposed  a  bearing  stiffness  matrix  which 
consists  of  kbrr,  kbr$  and  k^Q  coefficients.  Some  of  the  key  features  of  his  model  are 
also  summarized  in  Reference  [22].  This  model  is  in  fact  a  subset  of  our  [K]bm  as 
shown  in  Chapter  II  of  this  report  He  incorporated  his  bearing  model  in  a 
system  study  using  a  commercial  structural  synthesis  program  [38],  However, 
based  on  our  study  we  have  inferred  that  he  incorrectly  synthesized  the  system 
model  given  the  plate  experimental  modal  data,  shaft  finite  element  model  and 
analytical  bearing  model.  Moreover,  an  error  was  found  when  he  converted  kbre 
and  kb00  coefficients  to  "effective  stiffness  coefficients"  which  he  claimed  to 
couple  the  shaft  bending  motion  to  the  plate  out-of-plane  motion.  Also,  this 
method  excludes  the  bearing  rotational  degree  of  freedom,  which  from  our  study 
was  found  to  be  important. 

3.3  ASSUMPTIONS  AND  OBJECTIVES 

Linear  discrete  vibration  models  of  the  generic  system  shown  in  Figure  3.1  arc 
used  to  incorporate  [K]^  and  to  characterize  the  vibration  transmission  through  rolling 
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element  bearings.  The  stiffness  coefficients  of  [K]bm  are  evaluated  using  the  analytical 
expressions  presented  in  Chapter  II  of  this  report.  Effect  of  the  gyroscopic 
moment  on  the  shaft  dynamics  is  not  included.  Since  the  bearing  system  is  statically 
indeterminate,  the  direct  stiffness  formulation  technique  is  used  to  obtain  the  system 
governing  equations  as  opposed  to  the  flexibility  formulation.  The  governing  equations 
for  the  system  vibration  model  can  be  given  in  the  matrix  form  as 

[M]{q(t)}a+[C]{q(t)}a+[K]{q(t)}a  =  {f(t)}a  (3.1) 

where  [M],  [C]  and  [K]  are  the  system  mass,  damping  and  stiffness  matrices 
respectively,  and  { q(t) } a  and  {f(t)}a  are  defined  as  the  generalized  alternating 
displacement  and  applied  load  vectors  respectively.  Due  to  the  linearity  of  the  vibrating 
system,  mean  shaft  loads  { f } bm  a°d  preloads  do  not  directly  affect  the  dynamic 
response  of  the  rotating  system  and  hence  are  excluded  from  equation  (3.1).  However, 
{f)bm  and  bearing  preloads  are  assumed  to  be  constant  to  ensure  a  time-invariant  [K]bm 
matrix  which  depends  only  on  these  mean  loads  or  on  the  mean  deflection  operating 
points.  Accordingly,  only  the  alternating  shaft  loads  {f(t)}sa  in  Figure  3.1  which 
represent  typical  machine  excitation  due  to  the  kinematic  errors,  mass  unbalances  and 
torque  fluctuations  are  included  in  the  forced  vibration  problem.  The  energy  dissipation 
associated  with  the  rolling  element  bearings  is  assumed  to  be  an  energy  equivalent 
viscous  damping  matrix  [C]b  =  0  [K]bm  where  c  is  the  Rayleigh  damping  matrix 
proportionality  constant.  Dynamic  instabilities  due  to  the  oil  whirl  phenomenon  and 
asymmetry  of  rotating  elements  [11,12]  are  clearly  beyond  the  scope  of  this  study  and 
hence  are  not  considered  here. 
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The  specific  objectives  of  this  chapter  are  to:  (i)  incorporate  the  proposed  bearing 
matrix  [K]bm,  developed  in  Chapter  II  of  this  report,  in  the  linear  discrete  vibration 
model  of  the  rotating  mechanical  equipment  as  described  by  equation  (3.1)  using  both 
the  lumped  parameter  and  dynamic  finite  element  methods,  (ii)  evaluate  the  dynamic 
stability  of  the  proposed  bearing  system  model  using  the  Liapunov's  second  method, 
(iii)  calculate  eigensolution  and  forced  harmonic  responses,  and  predict  vibration 
transmission  through  rolling  element  bearings  for  three  example  cases,  (iv)  demonstrate 
the  advantages  of  our  formulation  over  the  existing  models  by  Kraus  et  al.  [24]  and 
White  [20],  and  (v)  validate  the  proposed  theory  by  comparing  analytical  prediction  with 
experimental  data  on  an  analogous  system. 

3.4  SYSTEM  GOVERNING  EQUATIONS 
3.4.1  Method  A:  Lumped  Parameter  Model 

The  proposed  bearing  matrix  [K]bm  can  be  easily  implemented  in  equation  (3.1). 
Note  that,  the  coupling  coefficients  of  [K]bm  provide  the  capability  to  predict  casing 
rigid  body  angular  0jca(t),  j=x,y,z,  and  translational  ujca(t)  motions  given  only  the 
unidirectional  transverse  shaft  forces.  Hence  we  can  couple  the  shaft  motions  to  the 
motions  of  a  casing  of  a  system  similar  to  Figure  3.1  but  with  rigid  shaft  and  rigid 
casing  using  a  lumped  parameter  model.  The  bearing  preloads  can  now  be  included  in 
the  mean  shaft  load  vector  {f)bm  by  a  direct  vector  addition  as  the  rigid  shaft  can  be 
assumed  to  be  a  single  lumped  mass  for  this  purpose.  An  alternating  displacement 
vector  [q(t)]a  =  {{q(t)  [q(t)}^  }T,  is  defined  where  (qW),*  =  {ujsa(t),0jsa(t)]T 

and  [q(t))ca  =  {ujca(t),0jca(t)J  ,  j=x,y,z,  are  the  shaft  and  casing  alternating 
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displacement  vectors  respectively.  The  governing  equations  of  motion  for  this  generic 
vibration  model  with  DOF  =  12  are  given  by  equation  (3.1)  with 
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where  the  stiffness  matrices  and  [K]v  pertain  to  the  bearing  and  mount 

respectively,  and  the  matrices  [M]s  and  [MJC  are  diagonal  shaft  and  casing  mass  matrices 
respectively;  each  matrix  of  dimension  6.  Specific  examples  of  this  method  along  with 
the  eigensolution  and  forced  response  studies  will  be  presented  in  Sections  3.7  and  3.8. 

3.4.2  Method  B:  Dynamic  Finite  Element  Formulation 

Consider  the  dynamic  finite  element  method  of  incorporating  [K]bm  in  equation 
(3.1)  especially  when  shaft  and  casing  plate  are  elastically  deformable  over  the 
frequency  range  of  interest.  This  method  is  different  from  the  lumped  parameter 
formulation  of  Section  3.4.1  which  assumes  non-compliant  shaft  and  casing.  For 
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example,  if  the  flexible  casing  plate  is  considered  to  be  very  large  compared  to  the 
bearing  dimensions,  then  the  bearing  nodal  point  on  the  shaft  can  be  coupled  to  only  one 
bearing  nodal  point  on  the  plate  as  shown  in  Figure  3.2a.  Accordingly,  the  present  form 
of  [K]bm  is  implemented  in  the  finite  element  model  as  a  generalized  stiffness  matrix  like 
the  lumped  parameter  model.  On  the  other  hand,  when  the  flexible  casing  plate 
dimensions  are  finite  and  of  the  order  of  bearing  dimensions,  then  several  bearing  nodal 
points  are  considered  as  shown  in  Figure  3.2b.  The  discretization  philosophy  here 
assumes  that  a  relative  displacement  vector,  given  by  the  difference  between  the 
averaged  displacement  vector  of  bearing  nodal  points  on  the  plate  and  the  displacement 
vector  of  a  bearing  nodal  point  on  the  shaft,  is  equivalent  to  the  actual  rigid  body  bearing 
motion.  Accordingly,  we  divide  the  bearing  stiffness  coefficients  equally  among  all  the 
generalized  stiffness  elements  connecting  the  bearing  nodal  points  on  the  plate  to  a  single 
bearing  nodal  point  on  the  shaft.  In  the  limit,  where  all  the  bearing  nodal  points  on  the 
plate  are  collapsed  to  a  single  nodal  point,  [K]bm  is  recovered  as  in  the  first  method. 
Our  finite  element  formulation  uses  conventional  structural  elements  typically  available 
in  commercial  software  programs  [39]  -  this  will  be  illustrated  in  Section  3.9.  Other 
features  of  this  method  are  similar  to  those  discussed  earlier  in  Section  3.4.1. 

3.4.3  Other  Methods 

Alternate  methods  of  incorporating  [K]bm  in  equation  (3.1)  such  as  finite 
difference  which  is  similar  to  method  B,  flexibility,  component  mode  synthesis  and 
transfer  matrix  formulations  are  also  possible.  In  the  flexibility  formulation,  the  bearing 
flexibility  matrix  can  be  obtained  by  inverting  a  subset  of  the  bearing  stiffness  matrix 
[K]bms  which  excludes  zeroes  corresponding  to  9Z  angular  direction  from  lK]bm.  In  the 
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Flexible  shaft 


(a)  Single  Bearing  Nodal  Point  on  the  Plate 


Figure  3.2  Discretization  method  for  implementation  of  [K]^  in  finite  element  model 
of  a  system  similar  to  Figure  3.1.  Here  {q(t)}sa  and  (q(i)}ca  are  the 
alternating  shaft  and  plate  displacement  vectors  at  bearing  location 
respectively.  Subscript  a  implies  alternating  component 
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transfer  matrix  method,  the  field  matrix  [T]b  for  a  bearing  can  be  easily  related  to 
[K]bms* 


[K] 

[I] 


-  1 
bms 


(3.3a,b) 


where  { {q(t)  }^,  (f(t)  }T  and  { {q(t)  }Ta,  (f(t)  )Ta  }T  are  now  the  state  vectors 
at  bearing  locations  on  the  shaft  and  casing  plate  respectively.  Equation  (3.3)  can  now 
be  integrated  with  transfer  matrices  of  the  shaft  and  plate  which  are  well  documented  in 
References  [40,41].  Direct  application  of  these  alternate  methods  [40-43]  are  beyond 
the  scope  of  this  paper  and  are  left  for  further  research. 


3.5  BEARING  SYSTEM  STABILITY 

The  stability  of  the  proposed  linear,  non-gyroscopic  model  of  a  bearing  system 
similar  to  Figure  3.1,  which  is  governed  by  equation  (3.1)  with  { f(t) } a  =  [0],  can  be 
determined  using  several  techniques  such  as  Liapunov's  stability  method,  Routh- 
Hurwitz  criteria,  or  from  the  direct  evaluation  of  system  eigenvalues.  Here  the 
Liapunov's  second  method  is  used  for  its  simplicity  when  applied  to  such  a  vibration 
model  [44,45].  If  the  system  matrices  [M],  [C]  and  [K]  of  equation  (3.1)  are  always 
symmetric  and  positive  definite,  then  the  system  is  asymptotically  stable  per  Liapunov. 
The  first  condition  is  directly  satisfied  since  [M],  [C]  and  [K]  are  symmetric.  Further, 
since  [M]  is  diagonal  and  consists  of  only  positive  entries,  it  is  clearly  positive  definite. 
For  [K],  the  positive  definite  test  can  be  performed  by  evaluating  its  principal  minor 


determinants  which  is  demonstrated  here  for  the  generic  lumped  parameter  model. 
Consider  the  decomposition  of  [K]  given  by  equation  (3.2)  into  a  product  of  3  matrices 


where  the  square  submatrix  [0]  of  the  appropriate  dimension  consists  of  only  zero 
entries,  [I]  is  an  identity  matrix  of  the  same  dimension,  and  [K]bm  and  [K]v  have  been 
defined  in  Section  3.4.1.  The  determinant  of  [K]  is  the  product  of  the  determinants  of 
the  three  matrices  on  the  right  hand  side  of  equation  (3.4) 

|k1  =  |[K]v||[K)J  (3.5) 

If  [  K]  I  >  0,  w j  =  1,2,... ,P  and  P  =  1,2,..., 12,  [K]  is  positive  definite.  Other 

wj  I 

principal  minor  determinant  with  Pel 2  can  be  obtained  by  excluding  the  stiffness 
coefficients  which  are  not  entries  in  the  principal  submatrix  of  equation  (3.5).  Since 
I  [  Kv]  |  =  kvx*kvy*kvz»kv0x»kYQy* kvQz >  0,  it  implies  that  equation  (3.5)  is  positive 

only  if  |  [  K 1  bm  |  is  positive.  We  may  recall  that  [K]^  has  zero  entries  in  the  last  row 

and  in  the  last  column  corresponding  to  torsional  9Z  angular  direction  which  forces  the 
bearing  system  to  be  semideftnite.  Now  define  a  new  matrix  [K]^  of  dimension  5  as 
a  subset  of  [K]bm  with  these  zeroes  excluded.  If  [K]bms  has  positive  principal  minor 
determinants,  then  this  system  is  dynamically  stable  because  it  will  consist  of  stable 
oscillations  superimposed  on  the  mean  shaft  rotational  motion  Qz  *  £2z(t).  Further,  it 
follows  that  [C]bms  which  is  proportional  to  [K]^  is  also  positive  definite  if  [K]bms  is 
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positive  definite.  The  resulting  equation  obtained  from  the  expansion  of  the  determinant 
of  [K]bms  in  terms  of  its  entries  kbwj,  wj=l,2,3,4,5,  is  given  as  follows  in  terms  of  the 
stability  functions  Oj,  j  =  1,2,3, 4. 


[K] 

bins 


=  2  04(55,34,45,35)  {kb24  04(1 1,23, 12,1 3)+kbl4  04(22, 13, 12,23)}+ 
2  <D4(44,35,45,34)  {kb25  04(ll,23,12,13)+kbl5  04(22, 13, 12,23) }  + 

2  04(24,35,25,34)  {kbl4  04(25, 13,23, 35)+kbl5  04(  12,34,24, 13)}  - 


kbi5  kb22  kj,^  <1>3(14,35,15,34)  +  kb24  kbi5  kbj4  03(23,35,33,25)  + 
2  kbl4  04(25,45,24,55)  04(13, 23,12,33)  +  {  kbl4  kb25  }2  kb33  + 


Oi(l,2)  02(45,34,35,44,55)  +  0^1,3)  02(45,24,25,44,55)  + 
kbl5  0^2,3)  03(14, 45, 15,44)  +kbll  {  04(24,35,25,34)  }2  + 
{kbi4}2  {  kb23  03(23, 55,35, 25)  }  + 

Oi(4,5)  {  kbll  0^2,3)  +  02(23, 13,12,33,22)  )  + 

2  kbl5  04(24,45,25,44)  04(23,13,12,33)  + 

2  kbl5  04(14,25, 15,24)  04(23,34,24,33)  (3.6a) 

Oi  (wj.wjj  )  =  kbWjWj  kbWjjWjj  —  { kbWjWjj }  (3.6b) 


02(wi,wn,win,wiv,wv)  —  2kbwj  kbwn  kbw^  kbWjV  {  kbWffl  }  kbWy  {  kbw^  }  (3.6c) 


^3  (wi>wii-wiii*wiv  )  =  2kbwi  kbWjj  -  kbWjjj  kbWjy  (3.6d) 

^4  (wi,wu,wm,wIV )  =  kbWj  kbWn  -  kbWm  kbWlv  (3.6e) 


where  wj,  j  =  I,II,III,IV,V,  are  the  dummy  variables  and  each  wj  may  represent  either  a 
single  number  or  a  set  of  two  numbers  in  equation  (3.6a).  The  principal  minor 
determinants  of  [K]bms  can  also  be  derived  from  equation  (3.6)  by  excluding  the 
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appropriate  stiffness  coefficients  which  are  not  entries  of  the  particular  principal 
submatrix.  Hence,  the  bearing  system  is  stable  if  each  principal  minor  determinants 
derived  from  equation  (3.6)  is  positive.  The  stability  of  the  proposed  bearing  model 
given  in  Chapter  II  of  this  report  can  now  be  verified  using  these  conditions. 
Inequalities  associated  with  the  stability  criteria  for  these  models  are  summarized  in 
Tables  3.1  and  3.2  for  ball  and  roller  bearings  respectively.  These  inequalities  arise  due 
to  the  fact  that  we  are  yet  to  impose  any  restrictions  on  these  stiffness  coefficients.  In 
Chapter  II  of  this  report,  these  stiffness  coefficients  are  given  as  functions  of 
bearing  kinematic  and  design  parameters,  and  hence  any  coefficient  can  not  assume  an 
arbitrary  value  as  it  is  related  uniquely  to  other  coefficients  through  these  parameters. 
Extensive  numerical  studies  performed  over  a  wide  range  of  these  parameters  have 
indicated  that  the  bearing  models  proposed  in  Chapter  II  are  indeed  stable  provided  the 
preloads  are  sufficiently  large  to  avoid  the  clearance  non-linearity.  Figures  3.3-3.5 
illustrate  examples  of  these  bearing  system  stability  studies  for  precision  rolling  element 
bearings  whose  design  data  are  given  in  Table  3.3.  In  all  of  these  figures,  the  stability 
functions  Oj  given  in  Tables  3.1  and  3.2  are  found  to  lie  within  the  stable  region. 

3.6  SYSTEM  RESPONSE 

The  eigensolution  of  the  linear,  non-gyroscopic  undamped  system,  formulated  by 
setting  {f(t)}a  =  (0)  and  [C]  =  [0]  in  equation  (3.1)  given  by  method  A  or  B,  yields  real 
valued  natural  frequencies  0)j,  j=l,2,3,...,  and  the  modal  matrix  [U]  =  [{<t>)i,  {<j>}2. 
•••»  •••]  for  the  stable  system.  Since  the  system  is  proportionally  damped,  the 

modal  damping  ratio  is  =  c  to? /2  and  the  damped  natural  frequency  is 
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Table  3. 1  Bearing  system  stability  criteria  for  the  proposed  ball  bearing  model 
(j  =  x  or  1,  y  or  2;  p  =  x  or  1,  y  or  2  but  p  *  j) 


Proposed  bearing  model^ 
see  Table  2.2 

Stability  criteria^ 

^•xx’^yy^zz’^OxOx’^flyOy’^zOp 

<t»j(3,p+3)  >  0 

kxx'^yy’kzz’kQx0x>kQy0y,kxgy,kyQx 

<X>,(1,5)  >  0  ;  <X>j(2,4)  >  0 

kxx>kyy>kZZ’k0x0x>k0y0y 

always  stable 

kxx,kyy»kzz>k0xOx’k0y0y,kpz 

Ol(p,3)  >  0 

kxx>kyy.kzz,kgx0x,k0y0y,kXy, 

k0x0y»kz0x,kz0y 

<D1(1,2)>0;01(3,4)>0; 
k33<l>1(4,5)-Ml>2(4  5,3  4,3  5,4  4,5  5)  >  0 

kxx^yy^zz^xOx’^yOy^xGx* 

kyOy>kjz’kz0p 

kllk22^i(3»4)+kj32kppk44+k142k22k33  >  0 ; 

O1(j,3)>0; 

{k330,(1.4)«1(2^)05(p  P+3J+3  j+3,p  p,4  4,5  5)} 
(kj3  k3  p+3)2 '  (®5(  j  j+3.P  P.l  1,2  2,j+3  j+3))  >  0 

^xx'kyy,kzz,k0x0x,k0y0y,kxy, 

^xz>^yz,^GxGy 

Oj(4,5)  >  0  ;  *,(1,2)  >  0  ; 
kn  G>1(1,2)+<I>2(2  3,1  2,1  3,2  2,3  3)  >  0 

kxx’kyy>kzz'k0x9x,k0y6y’kx0y> 

ky0x,kjz,kz0p 

^(Pd+S)  >  0 ;  <D1(  j,3)  >  0  ; 
kjj0!(3,p+3)+02(3  p+3j  3j  p+3,3  3, p+3  p+3)  >  0 

all  non-zero  except  8Z  terms 

I  [K]  |>0 

1  bmsl 

t  Here  the  subscript  b  which  implies  bearing  has  been  omitted  for  brevity, 
tt  Stability  functions  <J>j  and  <J>2  are  defined  by  equations  (3.6b)  and  (3.6c)  and  d>5  is 
given  by  <J>5  (  wj,w1I,wm,wIV,wv  )  =  {kbWj}2  kbwn  -  kbwm  kbWiy  kbWy 
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Table  3.2  Bearing  system  stability  criteria  for  the  proposed  roller  bearing  model 
(j  =  x  or  1,  y  or  2;  p  =  x  or  1,  y  or  2  but  p  *  j) 

Proposed  bearing  model^ 
see  Table  2.3 

Stability  criteria^ 

kii 

always  stable 

^xx’^yy»^zz’^0x0x,k0y0y,kx0y>^y0x 

d>1(l,5)>0;01(2,4)>0 

^zz’^OxOx’^OyOy 

always  stable 

^zz*^0x0x,k0y0y,kz0j 

®!(3j+3)  >  0 

lr  V  If 

Kxx»ILyy*Kxy 

Oj(U)>0 

^■xx'^yy’^zz’^GxQx'^OyOy’^xBy* 

<Dl(pj+3)  >  0  ;  <I>!(  j,3)  >  0  ; 

ky0x’kjz,kz0p 

kjjOj(3,p+3)-Hl>2(3  p+3j  3j  p+3,3  3,p+3  p+3)  > 

^zz'^0x0x»k0y0y,k0xOy,kz0x,kz0y 

0,(3,4)  >  0 ; 

k33<D1(4,5)+<I>2(4  5,3  4,3  5,4  4,5  5)  >  0 

all  non-zero  except  8Z  terms 

|[K]  |>0 

1  bmsl 

t  Here  the  subscript  b  which  implies  bearing  has  been  omitted  for  brevity, 
tt  Functions  <1^  and  $>2  are  defined  by  equations  (3.6b)  and  (3.6c) 


Stability  Function  (MN3/m)  Stability  Function  (MN2) 


STABLE 

S2(ball) 

Stroller) 

"" 

SI  (roller) 

sifham 

Axial  Deflection  (mm) 


Figure  3.4  Plot  of  stability  criteria  functions  for  ball  and  roller  bearings  subjected  to 
mean  bearing  axial  deflection  8^.  Here,  Si  =  d>i(2,4)  and  S4  =  G>j(l,5). 


Angular  Deflection  (rad) 
(b) 


Figure  3.5  Plot  of  stability  criteria  functions  for  ball  and  roller  bearings  subjected  to 
mean  angular  misalignment  Pym.  (a)  SI  =  <Dj(2,4)  and  S2  =  <!>j(l,3).  (b) 
S3  =kbll  <Pj(3,5)+<1>2(35,13,15,33,44) 
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Table  3.3  Design  parameters  for  typical  ball  and  roller  bearings  used  for  system  studies 


Parameters 

Ball  type 

Roller  type 

Load-deflection  exponent  n 

3/2 

10/9 

Load-deflection  constant  K„  (N/mn) 

1.0  E9 

1.0  E8 

Number  of  rolling  element  Z 

12 

14 

Radial  clearance  rL  (mm) 

0.00005 

0.00175 

Pitch  diameter  (mm) 

40.05 

38.00 

Ao  (mm)  t 

0.05 

- 

Unloaded  contact  angle  Oq 

40° 

15° 

t  Unloaded  distance  between  inner  and  outer  raceway  groove  curvature  centers 
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given  by  co  jd  =  co  - 1?.  .  Free  vibration  response  due  to  the  initial  conditions  is 

not  considered  as  only  the  steady-state  particular  solution  corresponding  to  sinusoidal  or 
periodic  load  vector  {f(t)}a  is  of  primary  interest  Define  excitation  by  the  Fourier  series 
expansion  as  {f(t)}a  =  X  { f } ap  eit0pt  where  cop  =  pco0,  o>0  is  the  fundamental 
frequency,  and  {f}ap  is  the  complex  Fourier  coefficient  load  vector.  The  steady-state 
particular  solution  {q(t)}a  is  given  by  the  normal  mode  expansion  technique  [40,42,43] 


as 


{q(oja«[mf5>}: 


(4>}T{  f) 

I  1  ap 


;  p=l,2...  (3.7a,b) 


An  alternate  approach  would  be  to  assume  the  harmonic  solution  for  the  alternating 

displacement  as  {q(t))a  =  X  {q}an  e1<0pt-  Substituting  this  and  {f(t)}a  definition  into 

P 

equation  (3.1),  we  get 

Adj[-  0)p[M]  +  X[K]1 

m*  ;  x=1+ia“p  <38a,b> 

where  the  operator  Adj  refers  to  the  adjoint  of  the  dynamic  stiffness  matrix.  Features  of 
this  method  are  summarized  in  References  [40,42,43].  Since  the  vibration  transmission 
across  the  bearing  is  the  primary  issue,  we  now  define  sinusoidal  load  transmissibility 
R(pco)  terms  between  two  arbitrary  locations  I  and  II  as 
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‘f  f  = 


wla^Vl 


WVT 


W,j  =  1,2,.. .,6 


(3.9) 


where  fwja  and  fjna  are  components  of  the  dynamic  load  vectors  at  two  arbitrary 
locations  I  and  II  respectively.  The  accelerance  A(oip)  and  mobility  V(<op)  transfer 
functions  with  motion  at  location  I  due  to  an  alternating  force  or  torque  fjIIa  applied  at 
location  II  on  the  shaft  are 


A 


)  = 


«wIa«V 


fjHa 


v(0)p)  =  A(WP)  ;  w,j  =  1,2 . 6 


(3.10a) 

(3.10b) 


where  q  wIa  is  a  component  of  the  acceleration  vector  at  location  I.  Other  frequency 
response  functions  can  also  be  defined  in  a  similar  manner  [40,42,43]. 

3.7  EXAMPLE  CASE  I:  RIGID  SHAFT  AND  PLATE  SYSTEM 
3.7.1  Vibration  Models 

Consider  the  mechanical  system  shown  in  Figure  3.6a  which  is  assumed  to  be 
freely  suspended  or  softly  mounted  such  that  [K]v  *  [0].  A  ball  bearings  (see  Table 
3.3)  with  constant  axial  preload  is  supporting  a  short  rigid  shaft  subjected  to  a  mean 
torque  Tzsm  *  Tzsm(t)  and  a  sinusoidal  radial  force  Frsa(t)  =  Frsal  ei<0ot  applied  very 
close  to  the  bearing.  A  lumped  parameter  model  with  DOF  =  12  is  proposed  in  Figure 
3.6b.  Conversely,  the  same  system  has  also  been  analyzed  by  Kraus  et  al.  [24]  using  a 
simple  vibration  model  with  DOF  =  1  as  shown  in  Figure  3.6c  with  only  k{,rr 


72 


(a) 


x 


Rigid  Casing 
Ball  Bearing 


Rigid  Shaft 


Free  boundary  condition 


Rigid  Shaft  { 

;-[q(t»sa  DOF  =  6 

(b) 

Ball  Bearing  { 

[C]b 

Rigid  Casing  { 

Mc 

H 

k{q(t)L  DOF =6 
r  “ 

J  f  xsa  © 

Rigid  Shaft  { 

ms 

]-{ 

u  xsa  CO 

(c) 

Ball  Bearing  { 

kbrr^  tjj 

Cbrr 

Rigid  Casing  { 

mc 

]H 

^  uxca(0 

Figure  3.6  Example  case  I:  freely  suspended  rigid  shaft,  ball  bearing  and  rigid  plate 
system  subjected  to  alternating  radial  force  F^t)  applied  at  the  shaft,  (a) 
Physical  system,  (b)  Proposed  multi-degree  of  freedom  vibration  model 
with  DOF  =  12.  (c)  Simple  model  by  Kraus  et  al.  [24]  with  DOF  =  1. 
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coefficient  The  bearing  stiffness  matrix  for  an  axially  preloaded  ball  bearing  in  Figure 
3.6b  has  non-negligible  stiffness  coefficients  kbxx,  kbyy,  kbzz,  k^^,  kb6yey,  kbx0y 
and  kby0x  which  are  functions  of  the  mean  axial  preload  as  given  in  Chapter  n  of  this 

report.  The  system  matrices  of  equation  (3.2)  can  be  modified  for  this  case  by 

suppressing  other  bearing  stiffness  coefficients  and  [K]v.  It  can  be  easily  observed 

from  equations  (3.1)  and  (3.2)  that  5  sets  of  uncoupled  differential  equations  exist  The 

simplest  three  sets  are  homogeneous  and  pertain  to  the  rigid  body  torsional  motions 

0^(0  of  the  shaft  and  casing,  and  axial  vibration  u2a(t)  of  the  shaft-casing  system  which 

are  of  no  interest  here.  The  remaining  two  sets  are  almost  identical  and  associated  with 

either  { uxa(t),0ya(t) ) T  or  {uya(t),0xa(t)}T  degrees  of  freedom  for  rigid  shaft  and  casing. 

If  the  coordinate  system  is  chosen  such  that  F^t)  line  of  action  coincides  with  the  x- 

axis,  then  the  steady-state  solution  to  the  set  of  differential  equations  in  terms  of 
T 

{ uya(t),0xa(t) )  is  trivial.  Hence,  the  problem  reduces  to  a  semi-definite  vibration 
system  with  DOF  =  4.  Accordingly,  rewrite  [M]s,  [M]c  and  [K]bm  in  equation  (3.2)  in 
terms  of  the  displacement  vector  {q(t)}a  =  {uxsaW’®  ysaW)  uxcaW>®ycaW  as 


(3.11) 


First  two  eigenvalues  corresponding  to  the  rigid  body  motions  in  x  and  0y  directions  are 
zero.  The  dimension  of  equation  (3.11)  is  further  reduced  to  DOF  =  2  by  defining 
relative  motions  8xa(t)  =  uxsa(t)  -  uxca(t)  and  Pya(t)  =0ysa(t)  -  0yca(t)  which  turn  out  to 
be  the  bearing  rigid  body  motions. 
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Ym  ms  ® 

0  Vs 


fS. 


xa 


yaj 


mc 

Ym  -  ms  +  mc 


(3.12a) 


(3.12b) 


It  may  be  noted  that  purely  translational  8xa(t)  model  by  Kraus  et  al.  [24]  shown  in 
Figure  3.6c  constitutes  a  subset  of  equation  (3.12)  with  kbxey  =  0.  Eigensolution  of 
equation  (3.12)  with  [C]=[0]  yields  the  following  natural  frequencies  tOj  and  modes 


B1±B2 


®i.  2  J  2ymm  sy{  Is  •  Bi=='lr'nmskbe  0 ^^i^bxx 


={  1  •  B3} 


2 

bx0. 


T 

{0}  « 


1,2 


B3  = 


Ymmskbeyey“YlIskbxx±B2 


2Yjk 


bx9. 


(3.13a,b) 


(3.13c) 

(3.14a) 


(3.14b) 


On  the  other  hand,  the  eigensolution  of  the  single  degree  of  freedom  system  is  given  by 
A  /  k.  A  T  f  1  ) 

t02  =  V  Y  m  ^  {0^  ~  |  /y  m  |  where  subscript  2  and  superscript  A  are 
chosen  to  indicate  that  this  solution  essentially  estimates  0)2  and  {<j>)2  in  equations 
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(3.13)  and  (3.14).  Since  o>2  does  not  include  kbxQy  and  kbQy0y,  a>2  <  as  evident 
from  Table  3.4  for  3  different  axial  preloads.  This  natural  mode  is  dominated  by  8xa(t) 
as  indicated  by  equation  (3.14).  And,  the  first  mode  {<}>}  i  which  is  predominantly 
pya(t)  is  also  affected  by  the  axial  preload. 


Table  3.4  Bearing  stiffness  coefficients  and  undamped  natural  frequencies  (Hz)  of 
example  case  I* 


Axial 

preload 

Fzbm(N) 

Bearing  stiffness  coefficients 

Proposed  model 
(DOF  =  2) 

Simple 

model 

(DOF=l) 

A 

to2 

kbxx  (N/m) 

kbx0y  (N) 

kbey0y(Nm) 

COj 

co2 

115 

1.84  E7 

-3.05  E6 

1.36  E4 

156 

372 

341 

190 

2.13  E7 

-3.12  E6 

1.70  E5 

191 

395 

367 

285 

2.43  E7 

-3.09  E6 

2.02  E5 

221 

416 

392 

t  Other  system  parameters  are:  ms=10.0kg,  IjsO.OlSkgm2,  Ym=0.4,ys=0.3,a=lE-6s. 


3.7.2  Bearing  Transmissibility 

The  forced  harmonic  response  of  equation  (3.12)  can  be  obtained  using  the 
dynamic  stiffness  approach  given  by  equation  (3.8) 
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^b0y6y^  -  's)  ^"xsale‘ 


5xa(i)  = - 

(kb0y0y^  “  WoVs)(kbxx^  “  woYmm  s)  -(kbxe/) 


(3.15a) 


P  (t)  = 

*ya  ' 


-kb 


(kb8y0yk  t°oYIIs)(kbxx^'  <0oYmms)  (kbx0y^) 


(3.15b) 


It  can  be  seen  from  equation  (3.15)  that  Fxsa(t)  not  only  excites  8^,(0  but  Pya(t)  as  well, 
which  is  not  predicted  by  Kraus  et  al.  [24],  The  steady-state  solution  for  this 
simple  model  is  given  by  the  following;  compare  it  with  equation  (3.15a). 


(3.16) 


Both  models  are  used  to  determine  the  load  transmissibility  magnitude  terms  R(co0) 
which  are  computed  using  {5xa(t),Pya(t)}  ,  [K]  and  [C].  Dynamic  bearing  force  Fxba(t) 
and  moment  Myba(t)  magnitudes  excited  by  the  shaft  force  Fxsa(t)  are  given  by  force 
transmissibility  RFxbt,Fxst(®o)  anc*  moment  transmissibility  ^Myba,Fxsa(®o) 

respectively. 


Rf  ,F  (“o) 

xba  xii 


^/(l  +  o2co2)2  B4  +  (1+  <t2(0q)  B 5 


J 


(l+a2co5)2B4+(l+a2(oJ)(B5  +  B6)+B7 


(3.17a) 


i&arrrii'flvf*' 
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RM  ,F 

yba’  xsa 


YmVs^o  J( l  +  <y2<»20)k 


2„2  x  ,.2 

bxG, 


A/(l  +  a2a)2)2B4+(l  +  o2CD20)(B5+B6)  +  B7 


B  =  ( k  k  -  k 

'  '  be  e  h*  bxey 


)' 


B5“t0o^I  Iskbxx(kbxxYIIsWo+  2kbxey  2kbxxkb0yey) 

B6-t0o^m  mskbeyey(k 


y”yv  b0yeyYmm  s“o+2kbx0y  2kbxxkbeyey 


) 


(3.17b) 

(3.17c) 

(3.17d) 

(3.17e) 


B7  -  2cooYmm  sY^s  2kbxxkb0y6y  t0o(YIIskbxx +Ymm  ^bOyOy)}  + 

ymm  s?i 1  s  { ym™  sYi 1  -  2kbx ey(1  -  02(0  o >}  (3>17f) 


Only  the  force  transmissibility,  as  given  below,  is  predicted  by  the  simple  model  [24]; 
compare  it  with  equation  (3.17a). 


(3.18) 


Figure  3.7  compares  equations  (3.17)  and  (3.18).  Our  model  predicts  higher 
R fxba*^xf ^  ^  l^e  simPle  model  due  to  the  additional  constraints  imposed 

by  [K]bm.  Also,  it  is  clear  that  the  simple  model  can  not  predict  dynamic  moment 
transfer  through  the  bearing.  The  bearing  transmissibility  functions  R(co0)  predicted  by 
our  model  for  3  different  axial  preloads  are  shown  in  Figure  3.8.  Note  that  the  resonant 


78 


Figure  3.7  Bearing  transmissibility  spectra  Rfwbt,Fxsa(wo)  for  example  case  I.  Here, 
Rl:  force  transmissibility  with  fwba  =  Fxba  and  R2:  moment  transmissibility 
with  fwba  =  Myba,  as  predicted  by  our  model  with  DOF  =  2  and  the  simple 
model  by  Kraus  et  al.  [24]  with  DOF  =  1. 


Moment  Transmissibility  (m)  Force  Transmissibility 
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Figure  3.8  Effect  of  mean  axial  preload  on  the  bearing  transmissibility  spectra 
^fwb*'Fxia(Q>°)  ^fined  by  equation  (3.17)  for  example  case  I.  (a)  Force 
transmissibility  with  f^  =  F^.  (b)  Moment  transmissibility  with  fwba  = 
Myba> 
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amplitudes  vary  depending  on  the  amount  of  axial  preload  and  the  resonant  frequencies 
increase  with  increasing  preloads  as  expected. 

3.8  EXAMPLE  CASE  II:  RIGID  SHAFT  AND  PLATE  SUPPORTED  ON 
FLEXIBLE  MOUNTS 

3.8.1  Vibration  Models 

The  physical  system  of  example  case  I  is  modified  to  include  flexible  mounts,  [K]v 
*  [0],  and  mean  radial  shaft  force  Frsm  *  Frsm(t)  and  F^ft)  =  F^e10^1  as  shown  in 
Figure  3.9a.  The  ball  bearing  (see  Table  3.3)  is  also  preloaded  in  the  axial  direction. 
This  is  modeled  using  lumped  parameter  theory  with  finite  mount  stiffness  coefficients 
kvj,  j  =  x,  y,  z,  0X,  0y  or  0Z  as  illustrated  in  Figure  3.9b.  Since  the  rigid  shaft 
assumption  still  holds,  Fism  =  F^  is  applied  directly  on  the  bearing  in  a  manner  similar 
to  the  axial  preload.  The  bearing  matrix  [K]^m  in  Figure  3.9b  has  non-negligible 
stiffness  coefficients  kbxx,  k^yy,  ^b0x0x*  ^b0y0y*  ^bx0y»  ^by0x  ^bxz  and  ^bz0x 
which  are  functions  of  the  mean  bearing  load  vector  (f  }bm~  (Asm  or  1116  mean  bearing 
displacement  vector  {q)bm  =  (q}Sm  *  {q}Cm  as  given  by  Chapter  II  of  this  report. 
Conversely,  White  120]  has  investigated  this  problem  using  a  simple  model  with  DOF  = 
2  as  shown  in  Figure  3.9c.  It  may  be  noted  that  his  bearing  system  model  did  not 
include  the  effect  of  axial  preload. 

Like  example  case  I,  the  governing  equations  (3.1)  and  (3.2)  can  be  modified  and 
reduced  to  4  uncoupled  sets  of  differential  equations.  The  first  two  sets  associated  with 
0ysa(t)  and  9yca(0  are  homogeneous.  The  third  set  is  similar  to  equation  (3.1 1)  but  with 
x  and  y  subscripts  interchanged,  and  two  mount  stiffness  coefficients  kvy  and  kv6x 

included.  However, it  is  still  homogeneous  and  therefore  has  only  trivial  steady-state 
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Figure  3.9  Example  case  II:  rigid  shaft,  ball  bearing  and  rigid  plate  system  supported 
by  flexible  mounts  and  subjected  to  alternating  radial  force  Fxsa(t)  applied  at 
the  shaft,  (a)  Physical  system,  (b)  Proposed  multi-degree  of  freedom 
vibration  model  with  DOF  =  12.  (c)  White's  vibration  model  [20]  with 
DOF  =  2. 
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solution.  The  final  set  corresponds  to  {q}a(t)  =  { uxsa(t),  u^tj,  0ysa(t),  uxca(t),  uzca(t), 
0yca(t)}T  with  [MJS,  [M]c,  [K]bni  and  [K]v  of  equation  (3.2)  reduced  to 


[K]  = 

bm 


[m.  °  01 

'm 

c  0 

0  ' 

[M]s  = 

0 

3 

to 

0 

9 

[Ml, 

0 

mc 

0 

Lo  0  i.j 

0 

0 

h. 

bxx 

k.  k  ‘ 

bxz  bx0  v 

"  k 

vx 

0 

0 

k. 

bxz 

1 

-O 

J 

M 

9 

IK)V  = 

0 

fc- 
>  < 

N 

0 

y 

0 

0 

k  a 

k 

V0 

bx9y 

bz0y  b0y0y 

'y  J 


(3.19a,b) 


(3.19c,d) 


The  vibration  model  by  White  [20]  may  be  formulated  by  retaining  only  two  equations 
corresponding  to  uxsa(t)  and  uxca(t)  and  excluding  all  bearing  and  mount  stiffness 
coefficients  except  for  kbxx  and  kvx  . 

Analytical  eigensolution  of  the  undamped  system  is  not  possible  since  it  requires 
solving  for  the  zeroes  of  a  6-th  order  polynomial  in  0)2.  Therefore,  this  problem  is 
solved  numerically  using  an  eigenvalue  routine  [42].  Using  the  same  system  parameters 
as  in  example  case  I  with  3  different  Fxsrn,  natural  frequencies  and  modes  are  found  as 
given  in  Tables  3.5b  and  3.5c  for  both  our  and  White's  models.  Corresponding  bearing 
mean  loads  (Fxbm>Fzbm>Mybm }  and  relevant  bearing  stiffness  coefficients  computed 
using  the  method  derived  in  Chapter  II  are  listed  in  Table  3.5a.  Tables  3.5b  and  3.5c 
indicate  that  only  the  first  and  fifth  modes  of  our  model  are  predicted  by  White's  model; 
here  superscript A  is  again  used  to  denote  estimation  based  on  the  simple  model.  The 
first  natural  frequency  tOj  predictions  by  both  models  are  very  similar.  But  (1)5 
prediction,  whose  mode  is  similar  to  the  second  mode  of  example  case  1,  indicates  a  few 
discrepancies.  White's  model  also  underestimates  this  natural  frequency  due  to  the 
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Table  3.5  Results  of  example  case  II  * 


(a)  Computed  bearing  mean  loads  and  stiffness  coefficients 


Shaft 
mean  load 

Fxsm 

(N) 

Bearing  loads  ^ 

Bearing  stiffness  coefficients 

Fzbm 

(N) 

|1 

^bxx 

(N/m) 

^)XZ  ^bx0y 

(N/m)  (N) 

^bzz  ^bz8y 

(N/m)  (N) 

^yGy 

(Nm) 

44 

198 

2.26E7 

5.84E7  -2.90E5 

8.50E7  -3.10E3 

1.71E4 

94 

216 

1.16 

2.64E7 

1.14E7  -2.33E5 

8.55E7  -1.36E4 

1.72E4 

122 

229 

1.38 

2.93E7 

1.40E7  -1.98E5 

8.60E7  -2.19E4 

1.74E4 

(b)  Undamped  natural  frequencies  (Hz) 


Mean 
shaft  load 

Fxbm(N) 

Proposed  model  (DOF=6) 

Simple  model  (DOF=2) 

®1 

©2 

<u3 

<04 

0J5 

“6 

A 

*n 

<3 

44 

93 

100 

111 

287 

350 

607 

97 

321 

94 

96 

100 

115 

300 

351 

614 

97 

345 

122 

97 

100 

117 

308 

364 

620 

98 

362 

(b)  Modes  of  vibration  ' 


Proposed  model  (DOF=6) 


Simple  model  (DOF=2) 


Wl 

W2 

W3 

We 

(ill 

{ i}5 

0.217 

0.031 

RVPSV 

mtm 

0.044 

0.218 

0.229 

0.203 

wmm 

BiVitfll 

0.236 

1.648 

0.084 

5.966 

-1.245 

-0.068 

0.171 

0.084 

-0.071 

0.107 

-0.031 

0.187 

-0.178 

-0.027 

0.194 

0.008 

-0.018 

-0.165 

0.167 

0.009 

0.917 

4.648 

3.295 

0.073 

t  Other  system  parameters  are:  ms=10.0  kg,  ls=0.025  kgm2  mc=15.0  kg, 
Ic=0.03  kgm2,  cr=lE-6  s,  kvx=lE7  N/m,  kvz=lE7  N/m,  kx6y=l  E5  Nm. 
ft  ^xbm  =  Fxsm 

ttt  These  are  for  mean  shaft  load  Fxsm  =  FXbm  =  94  N. 
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incomplete  bearing  stiffness  model  employed.  Here  again  our  predicted  modes  include 
8xa,  5za  and  pya  displacements  of  shaft  and  casing  components  which  are  not 
considered  by  White’s  model. 

3.8.2  Frequency  Response 

The  forced  harmonic  response  solution  is  also  determined  numerically  using  the 

dynamic  stiffness  approach  outlined  in  Section  3.6.  Driving  point  (with  j  =  s)  and  cross 
point  (with  j  =  c)  accelerance  spectra  i5‘qwjaFxs#(£0o) 316  given  in  Figure  3.10  with  qwja  = 

uxja>  uzja  or  9yja  f°r  both  models.  Here,  we  observe  that  White's  [20]  model 
overestimates  the  magnitudes  of  the  accelerance  and  can  not  predict,  unlike  our  model, 
axial  u^ft)  and  0ya(t)  angular  motions  on  the  shaft  and  casing.  Figure  3.1 1  shows  the 
bearing  transmissibility  spectra  Rfwba,Fxsa(W0)  for  fwba  =  Fxba,  Faa  or  Myba  which 

indicate  that  transmissibilities  corresponding  to  F^  and  Myba  are  not  predicted  using 
White's  model.  Such  loads  also  serve  as  mechanisms  for  vibration  transmission 
through  the  bearing  to  the  casing,  in  addition  to  Fxba.  The  mount  transmissibility 
spectra  Rfwva,Fxsa(wo)  as  shown  in  Figure  3.11  indicate  that  Fzva  and  Myva  are  also 

transmitted  to  the  mounts  in  addition  to  Fxva  due  to  casing  motions  in  x,  z  and  0y 
directions.  The  effect  of  mean  radial  bearing  force  Fxbm  on  the  load  transmissibilities 
R(co0)  through  the  bearing  is  shown  in  Figure  3.12  for  3  different  mean  loads.  We 
again  observe  that  the  resonant  amplitudes  and  frequencies  are  mean  load  dependent 
through  [K]bm.  Similar  trends  show  that  an  increase  in  mean  load  raises  the  resonant 
frequencies  although  the  effects  are  not  as  pronounced  as  those  found  in  example  case  I. 


Magnitude  (I/kg)  Magnitude  (lAg) 
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Figure  3.10  Accelerance  spectra  Aqw  a a((°b)  f°r  example  case  II  as  predicted  by  our 

formulation  and  White's  model  [20].  (a)  Driving  point  accelerance  with 
9wja  =  uxsa-  0>)  Cross  point  accelerance  with  qwja  =  uxca.  (c)  Accelerance 
with  =  uzsa  for  shaft  and  qwja  =  u^a  for  casing,  (d)  Accelerance  with 
qwja  =  Oysa  for  shaft  and  qwja  =  0yCa  for  casing. 


Force  Transmissibility  Force  Transmissibility 
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Figure  3.1 1  Bearing  Rfwbt,FXM(ft)o)  and  mount  Rf^.F^/^o)  transmissibility  spectra 

for  example  case  II  as  predicted  by  our  formulation  and  White's  model 
[20].  (a)  Bearing  force  transmissibility  with  fwba  =  Fxba.  (b)  Mount  force 
transmissibility  with  fwva  =  Fxva.  (c)  Force  transmissibility  with  fwba  = 
FZba  for  bearing  and  fwva  =  Fzva  for  mount,  (d)  Moment  transmissibility 
with  fwba  =  Myba  for  bearing  and  fwva  =  Myya  for  mount 


Moment  Transmissibility  (m)  Force  Transmissibility  Force  Transmissibility 
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Figure  3.12  Effect  of  the  mean  radial  bearing  force  Fxbm  on  the  bearing  transmissibility 
spectra  Rfwbi,FXfa(<Bb)  for  example  case  n.  (a)  Force  transmissibility  with 
fwba  =  FXba-  (b)  Force  transmissibility  with  fWba  =  Fzba.  (c)  Moment 
transmissibility  with  fWb8  =  Mybf, . 
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3.9  EXAMPLE  CASE  HI:  EXPERIMENTAL  STUDY 

3.9.1  Physical  setup 

The  final  example  case  examines  the  experimental  setup  of  Lin  [17]  as  shown  in 
Figure  3.13a.  This  system  is  similar  to  Figure  3.1  and  consists  of  a  159mm  long  x 
25mm  diameter  non-rotating  shaft  supported  by  two  rolling  element  bearings  of  25mm 
bore  x  51mm  outer  diameter.  One  is  supported  on  a  rectangular  plate  of  approximate 
dimensions  762mm  x  457mm  x  9mm  and  the  second  is  rigidly  connected  to  the  base. 
The  plate  is  also  bolted  to  a  massive  base  structure.  Excitation  force  Fys  which  consists 
of  a  mean  Fysm  =  445N  via  a  preloaded  spring  and  an  alternating  Fysa(t)  component 
applied  transversely  at  the  free  end  of  the  shaft  using  a  vibrating  shaker.  Driving  and 
cross  point  accelerance  spectra  are  measured  at  the  shaft  and  on  the  plate  respectively. 
Further  details  of  this  experiment  are  summarized  in  Reference  [17]. 

3.9.2  Bearing  Analysis 

Initially,  only  the  static  analysis  is  performed  to  obtain  [K]^  for  this  experimental 

system  using  the  method  proposed  in  Chapter  II  of  this  report.  The  static  analysis 

neglects  plate  flexibility;  this  assumption  is  valid  since  the  bearing  mean  loads  are 

sufficiently  low  and  do  not  deflect  the  plate.  The  shaft-bearing  system  is  statically 

indeterminate  as  shown  in  Figure  3.13b.  The  mean  force  Fysm  on  the  shaft  produces 

T 

mean  bearing  load  vector  {f)bm  =  (0-Fybm>Fzbm>Mxbm»C)}  which  depends  on  the 
mean  bearing  displacement  vector  {q}bm  =  {0,8ym,8zm,Pxm,0}  .  The  proposed 
bearing  matrix  [K]bm  includes  stiffness  coefficients  kbxx,  kbyy,  kbzz,  kb0x0x,  kb0y0y, 
kbx0y,  kby0x,  kbyz  and  kbz0x  which  are  direct  functions  of  {q)bm.  In  contrast,  the 


y 
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figure  3.13  Example  case  III.  (a)  Schematic  of  the  experimental  setup  consisting  of  an 
overhung  shaft,  2  ball  bearings  and  a  rectangular  plate  [  17].  (b)  Static 
model  of  shaft  and  bearings  used  for  computing  [K]^.  (c)  Finite  element 
model  of  the  experimental  system  shown  in  (a).  The  generalized  stiffness 
matrix  elements  for  the  bearing  model  are  shown  by  dashed  lines, 
connecting  1  node  on  the  shaft  to  4  nodes  on  the  plate. 
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conventional  bearing  models  include  only  kbxx,  kbyy  and  k^,T  coefficients.  From  Figure 
3. 1 3b,  the  force  and  moment  equilibrium  equations  for  this  system  are 

Fybjm  +  Fyb2m  -  ^ysm  =  ®  (3.20a) 

^ybim  +  ^yb2m  —  ^yb2m^I  +  Eysm^l  +  ^2)  =  ®  (3.20b) 

^zbim  +  ^zb2m  =  ®  (3.20c) 

Since  the  shaft-bearing  system  is  statically  indeterminate,  bending  theory  for  the  shaft 
and  rigid  body  motion  constraint  in  the  z  direction  are  used  to  estimate  stiffness 
coefficients 


El  Sy2m  —  El  Sylin  +  Fybim^lV^)  +  Mybjm^iV^) —  0 

(3.21a) 

El  pX2m  -  El  PXim  +  Fybl„ +  Myb^/j  =  0 

(3.21b) 

=  0 

(3.21c) 

Additionally,  6  nonlinear  algebraic  equations  defined  by  the  mean  bearing  load- 
displacement  relations  as  given  in  Chapter  II  are  required.  These  nonlinear  algebraic 
equations  are  solved  using  Newton-Raphson  method  [36,37].  Since  mean  loads  on 
each  bearing  are  sufficiently  large,  the  bearing  stiffness  coefficients  for  both  bearings  are 
almost  identical,  as  listed  in  Table  3.6  along  with  other  system  parameters. 
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Table  3.6  Design  and  estimated  parameters  for  two  identical  rolling  element  bearings 
used  in  example  case  Ill 


Load-deflection  exponent  n=3/2 

Ao^=0.05  mm 

kbzz=1.72E8  N/m 

Loaddeflection  constant  K„=6.92E9  N/mn 

/j=41  mm 

kbx0y=-2.56E5  N 

Number  of  rolling  element  Z=10 

/2=84mm 

kby0x=3.52E5  N 

Radial  clearance  rL=5.0E-5  mm 

kbxx  =1.44E8  N/m 

kbZQxi=4  02E5  N 

Pitch  diameter=38.1  mm 

kbyy  =3.69E8  N/m 

kbexex=4.19E4  Nm 

Unloaded  contact  angle  ao=0° 

kbyz=2.04E8  N/m 

kb0y0y=l.O2E4  Nm 

t  Unloaded  distance  between  inner  and  outer  raceway  groove  curvature  centers 


3.9.3  System  Study 

We  incorporate  the  proposed  rolling  element  bearing  stiffness  matrix  [K]bm  in  a 
finite  element  model  which  includes  shaft  and  plate  dynamics,  using  the  formulation 
given  in  Section  3.4.2.  The  finite  element  model  shown  in  Figure  3.13c  is  implemented 
with  a  commercial  software  [39].  The  shaft  component  is  modeled  using  2  noded 
Timoshenko  beam  elements  with  axial  degrees  of  freedom  in  addition  to  the  bending 
motion.  The  plate  model  is  constructed  using  4  noded  quadrilateral  plate  elements  with 
shear  deformation  and  rotary  inertia  effects.  Each  node  has  3  translational  and  3 
rotational  degrees  of  freedom.  Four  generalized  stiffness  matrices  corresponding  to  the 
first  bearing,  each  matrix  being  equivalent  to  1/4  [K]bm,  are  used  to  couple  the  single 
shaft  node  to  4  plate  nodes.  The  second  bearing  connects  one  end  of  the  shaft  to  a 
grounded  node.  The  boundary  conditions  for  the  plate  along  the  perimeter  are  chosen  to 
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be  a  combination  of  ideal  clamps,  Uxa(t)  =  uya(t)  =  ua(t)  =  G^ft)  =  0ya(t)  =  9za(t)  =  0, 
and  simple- supports  uza(t)  =  0  as  shown  in  Figure  3.13c  in  order  to  represent  the 
physical  model  as  much  as  possible.  Here,  the  energy  dissipation  is  assumed  to  be 
given  by  the  modal  damping  ratio  £  =  0.03.  A  sinusoidal  force  Fysa(t)  =  Fysaie^is 
applied  at  one  end  of  the  shaft  to  simulate  the  experiment 

Over  the  frequency  range  of  400Hz  to  2000Hz,  Figure  3.14  compares  the  driving 
point  accelerance  spectra  AUysatFysa(<o0)  =  u  ysa/Fysa.  The  simple  theory  shown  here 

represents  the  conventional  way  of  modeling  bearings  while  the  other  features  are 
exactly  the  same  as  in  proposed  model.  Our  predictions  match  measured  spectra  very 
well.  Conversely,  the  simple  model  predicts  slightly  higher  accelerance  amplitude  and 

lower  resonant  frequency  in  the  vicinity  of  800Hz  due  to  the  incomplete  bearing  model 
used.  Cross  point  accelerance  spectra  AUzcgiFysji(0)o)  =  u  zca/Fysa  are  shown  in  Figures 

3. 15a  and  3. 15b  where  u  is  measured  for  2  different  locations  on  the  plate  as  shown 
in  Figure  3.13c,  and  excitation  Fysa(t)  is  once  again  applied  transversely  at  the  shaft. 
Here,  each  predicted  accelerance  spectrum  has  been  averaged  over  4  points  in  the 
immediate  vicinity  of  the  measured  location.  Reasonable  comparisons  between  the 
proposed  model  and  experiment  are  seen.  Here,  the  discrepancies  are  primarily  due  to 
physical  setup  complexities  and  the  limitations  associated  with  the  finite  element  model 
in  describing  some  of  these.  In  Figure  3.15,  the  simple  model  is  not  included  because  it 
predicts  exactly  zero  out-of-plane  or  flexural  motion  of  the  plate.  Next,  the  cross  point 
mobility  level  Ly  is  defined  by  averaging  mean  square  mobility  spatially  over  the  entire 
plate  and  over  a  frequency  bandwidth  Aco.  This  level  is  directly  related  to  the  structure- 
borne  noise  or  vibratory  energy  transmitted  through  the  bearing. 


Magnitude  (1/kg) 
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Figure  3.14 


Driving  point  accelerance  spectra  Au  j;  (o0)  yielded  by  the  proposed 
model,  simple  model  and  experiment  l?y  Lm  [17]  for  example  case  HI. 
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•*  M 
u  u 
zca  zca 

F  F* 
xsa  xsa/J 

dB  re  cV2^  =  1 .0m2/N2s2  (3.22) 

where  ( )*  implies  the  complex  conjugate,  Sc  is  the  plate  surface  area  and  Re{  }  implies 
the  real  part  of  the  complex  number.  Table  3.7  compares  Ly  predictions  by  proposed 
and  simple  models  with  experimental  data.  It  can  be  seen  from  this  table  that  the 
proposed  model  predicts  the  experiment  quite  well  and  the  simple  model  fails  to  predict 
any  plate  vibration.  We  can  therefore  conclude  that  our  model  is  indeed  valid  for 
vibration  transmission  analyses. 


Lv=10I°g10{<v2zct>}  =  l0I°g 


10 


2  S.Ato 

c  Aco  S 


Table  3.7  Predicted  and  measured  cross  point  mobility  level  Ly  as  defined  by  equation 
(3.22) 


1/3  Octave  band 
center  frequency  (Hz) 

Experiment  (dB) 
[Lin] 

Proposed 
model  (dB) 

Simple 
model  (dB) 

400 

-102 

-105 

-OO 

500 

-92 

-96 

-oo 

630 

-95 

-94 

•OO 

800 

-88 

-97 

•OO 

1000 

-87 

-95 

-OO 

1250 

-97 

-108 

-OO 

1600 

-108 

-115 

-OO 

2000 

-106 

-107 

-OO 
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3.10  CONCLUDING  REMARKS 

A  new  mathematical  model  for  precision  rolling  element  bearing  has  been 
developed  and  incorporated  in  linear  system  dynamic  models  using  lumped  parameter 
and  finite  element  modeling  techniques  for  the  vibration  transmission  studies  of  a 
generic  single  shaft-bearing-plate  system.  Stability  studies  indicate  that  the  bearing 
system  is  dynamically  stable  for  most  of  the  practical  designs.  Through  3  example  cases 
including  one  experimental  study,  we  have  shown  that  our  proposed  vibration  model  is 
clearly  superior  to  the  models  currently  available  in  the  literature.  The  current  models 
tend  to  underestimate  the  resonant  frequencies  and  force/moment  transmissibilities,  and 
overestimate  the  accelerance  amplitudes  as  compared  to  our  models.  The  proposed 
model  also  predicts  how  the  vibratory  bending  motion  on  the  shaft  is  transmitted  to  the 
casing,  illustrated  through  coupling  coefficients  of  the  proposed  bearing  stiffness  matrix 
[K]bm.  Finally,  the  forced  response  trends  indicated  that  increase  in  the  mean  bearing 
loads  increases  system  resonant  frequencies.  We  are  extending  this  model  to  predict 
vibration  transmission  in  rotating  equipment  with  multiple  shafts,  bearings  and  gears. 
Other  applications  are  evident  as  our  theory  is  general  in  nature.  However,  it  is 
restricted  to  linear  systems.  Bearing  non-linearities  are  being  examined  in  a  parallel 
study  [46]. 


CHAPTER  IV 


GEARED  ROTOR  SYSTEM  STUDIES 

4.1  INTRODUCTION 

The  focus  of  this  chapter  is  on  the  dynamic  analysis  of  geared  rotor  system  shown 
in  Figure  4.1  which  includes  a  spur  gear  pair,  shafts,  rolling  element  bearings,  motor, 
load,  casing  and  flexible  or  rigid  mounts.  For  this  purpose,  discrete  vibration  models 
are  developed  and  used  to  predict  vibration  transmission  through  bearings.  Also,  the 
effects  of  bearing,  casing  and  mount  dynamics  on  the  internal  rotating  system  dynamic 
characteristics  are  investigated.  Traditional  analysis  approaches  [15,25-29]  in  the  gear 
dynamics  area  concentrate  on  the  internal  rotating  system  and  exclude  dynamic  effects  of 
casing  and  flexible  mounts.  Moreover,  simple  bearing  models  are  typically  used  which 
assume  either  ideal  boundary  conditions  on  the  shaft  or  translational  spring  elements. 
Some  of  the  limitations  associated  with  current  bearing  models  have  been  discussed 
thoroughly  in  Chapters  II  and  III  of  this  report.  For  instance,  simple  bearing 
formulations  can  not  explain  how  the  vibratory  motion  may  be  transmitted  from  the 
rotating  shafts  to  casing  and  support  structures,  and  moreover  can  not  predict  the  effects 
of  bearing,  casing  and  mount  dynamics  on  the  internal  rotating  system  adequately. 
Chapter  II  of  this  report  also  presents  a  new  mathematical  model  for  the  precision 
rolling  element  bearings  to  clarify  this  issue  qualitatively  and  quantitatively,  and  Chapter 
III  utilizes  the  proposed  bearing  stiffness  matrix  [K]^  formulation  to  analyze  the 
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Rigid  or  Flexible  Casing 


Figure  4. 1  Schematic  of  a  geared  rotor  system  with  flexible  shafts  and  rigid  spur  gear 
and  pinion,  and  supported  by  4  rolling  element  bearings  on  a  flexibly 
mounted  casing.  The  shaft  at  the  motor  end  is  driven  at  mean  speed 
which  in  turns  drives  the  load  at  mean  speed  fl  zL<  This  physical  system  is 
symmetric  about  a  plane  intersecting  both  the  driving  and  driven  shafts. 
The  y  axis  is  assumed  parallel  to  the  gear  mesh  line  of  action  in  the  pressure 
angle  direction  at  mesh  point.  This  schematic  is  used  for  example  case  II 
and  III. 
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vibration  transmission  problem  in  a  generic  shaft-bearing-casing-mount  system.  Our 
bearing  formulation  has  already  been  validated  through  several  example  cases. 

This  study  extends  the  proposed  bearing  formulation  of  Chapters  II  and  III  to 
examine  vibration  transmissibility  in  rotating  mechanical  equipment  through  two  generic 
systems  which  are  a  spur  gear  pair  drive  and  a  single-stage  rotor  system  with  rolling 
element  bearings  as  shown  in  Figures  4.1  and  4.2  respectively.  It  should  be  noted  that 
Figure  4.2  can  be  treated  as  a  special  case  of  Figure  4.1-  this  will  be  illustrated  later. 
The  physical  systems  are  assumed  to  be  driven  at  mean  speed  £2  2  =  £2^  in  Figure 

4. 1  and  £2  in  Figure  4.2;  subscript  M  refers  to  the  motor  end.  Each  shaft  is  supported 

by  two  identical  rolling  element  bearings  which  are  modeled  as  stiffness  matrices  [K]bm 

of  dimension  6  as  proposed  in  Chapter  II  of  this  report.  These  generic  systems  may 

be  excited  by  motor  and/or  load  torque  fluctuations,  rotating  mass  unbalances  or  gear 

kinematic  transmission  error  e(t)  [29,47];  here  e(t)  is  defined  as  the  deviation  of  the 

relative  gear-pinion  angular  position  from  its  relative  ideal  location.  Typical  excitation  at 

T 

the  shaft  is  defined  by  an  alternating  load  vector  { f(t) } §3  =  {Fj^O/Tj^t)}  =  {f(t)}s  - 
{ 1  )smi  j=x,y,z,  where  Fj^t)  and  Tj^t)  are  the  alternating  force  and  torque  respectively, 
{f(t)]s  is  the  total  load  vector  of  dimension  6,  { f }  sn,  represents  the  mean  load  vector, 
and  superscript  T  implies  the  transpose.  In  the  case  of  the  geared  drive,  (fft)^  consists 
of  e(t)  at  the  mesh  point.  Additionally,  in  the  case  of  the  single-stage  rotor  system, 
transverse  forces  F^t)  and  Fysa(t)  due  to  mass  unbalances  are  of  interest.  Note  that  in 
the  dynamic  analysis  {f}sm  and  bearing  preloads  are  not  included  as  they  do  not  appear 
in  the  governing  equations  of  the  linear  vibration  model  but  are  used  for  computing 
[K]^.  Other  effects  such  as  bearing  coupling  coefficients,  motor  and  load  inertia,  and 


Figure  4.2  Schematic  of  the  single  stage  rotor  system  with  flexible  shaft  rotating  at 
mean  speed  Qz  and  supported  by  2  identical  rolling  element  bearings  on  a 
flexibly  mounted  rigid  casing.  This  physical  system  is  symmetric  about  a 
plane  intersecting  the  shaft  axial  axis  and  parallel  to  the  x-z  plane.  This 
schematic  is  used  for  example  case  I. 


casing  and  mount  dynamics  are  considered  as  the  emphasis  is  on  the  overall  system 
behavior. 

4.2  LITERATURE  REVIEW 

The  available  literature  on  bearing  models  commonly  used  for  internal  rotor  and 
gear  dynamic  analyses  has  already  been  discussed  in  Chapters  II  and  III  of  this 
report.  Although  there  are  numerous  publications  [10-13,15,25-29]  on  the  modal 
analyses  of  geared  drives  and  single  or  multiple-stage  rotor  systems,  very  little  has  been 
reported  directly  on  the  force  transmissibility  through  bearings,  and  the  dynamic  effects 
of  bearing,  casing  and  mounts  on  the  internal  rotating  system. 

4.2.1  Casing  and  Mount  Dynamics 

A  comprehensive  review  of  the  available  vibration  studies  of  casing  and  mounts, 
without  the  internal  rotating  system,  has  been  given  by  Lim  and  Singh  in  1987  [6]. 
Other  related  publication  not  included  in  Reference  [6]  are  by  Gaul  and  Mahrenholtz  [48] 
in  1984,  and  Smith  [49]  in  1988  who  have  developed  lumped  parameter  models  of 
flexible  machine  foundations,  excluding  the  internal  rotating  system  dynamics.  Their 
studies  report  rigid  body  translational  and  rotational  vibration  responses  of  the  casing 
when  excited  by  forces  and/or  moments.  However,  these  studies  are  not  specifically 
directed  towards  geared  rotor  systems,  and  dynamic  interactions  between  internal 
rotating  system  and  casing-mount  system  are  not  incorporated.  In  addition,  Lim  et  al. 
[50]  in  1989  have  also  performed  a  dynamic  finite  element  analysis  of  a  helicopter 
gearbox  excluding  the  internal  rotating  system.  Results  again  show  that  casing  rigid 
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body  vibration  modes  are  dominant  over  the  lower  frequencies  for  a  flexibly  mounted 
casing,  but  numerous  casing  plate  elastic  modes  are  observed  at  higher  frequencies. 

4.2.2  Gear  Dynamics 

Current  gear  dynamic  models  include  only  the  internal  rotating  system  and  simple 
bearing  models,  and  typically  exclude  casing  and  mount  flexibilities  [15,25-29].  In  few 
instances  casing  and  mounts  have  also  been  included  ,  but  restricted  to  cases  where  only 
purely  radial  and/or  axial  force  on  the  bearing,  and  unidirectional  vertical/horizontal  rigid 
body  motion  on  the  casing  are  modeled,  as  summarized  in  Reference  [6].  Such  models 
still  do  not  explain  how  the  vibration  is  transmitted  from  the  gear  mesh  to 
the  casing  and  into  the  machine  foundation  as  witnessed  in  previous  experiments 
[9,16,171. 

4.2.3  Rotor  Dynamics 

Similar  to  the  geared  rotor  models,  most  of  the  existing  rotor  dynamic  models 
concentrate  on  the  internal  rotating  system  and  address  issues  related  to  stability,  critical 
speeds  and  rotating  mass  unbalance  response  [10-13],  But  a  few  investigators  have 
included  support  flexibility  models  [51-53J.  For  instance  Lund  and  Wang  [51]  in  1986 
proposed  an  impedance  matching  approach  to  reduce  the  large  degrees  of  freedom 
(DOF)  required  in  such  models.  Using  one  example  case,  they  reported  that  support 
foundation  has  little  influence  on  the  internal  rotor  resonances,  but  then  cautioned  that 
this  might  not  be  true  for  other  systems.  Vance  et  al.  [52]  in  1987  incorporated 
measured  support  foundation  parameters  in  a  transfer  matrix  model  of  the  rotor  dynamic 
system,  and  concluded  that  the  omission  of  support  flexibility  may  miss  some  of  the 


rotor  critical  speeds  of  interest.  Earles  et  al.  [53]  in  1988  reported,  based  on  a  finite 
element  model  of  the  support  foundation,  that  the  rotor  critical  speeds  shift  by  including 
support  flexibility.  In  all  of  the  above  mentioned  studies,  the  overall  dynamic  behavior 
of  the  rotor  system  is  only  investigated  partially.  Moreover,  only  a  subset  of  journal 
bearing  [K]bm  coefficients  is  included  which  can  not  predict  completely  the  vibration 
transmission  through  bearings. 

4.3  ASSUMPTIONS  AND  OBJECTIVES 

Linear  discrete  vibration  models  of  the  geared  and  single-stage  rotor  systems 
shown  in  Figures  4.1  and  4.2  are  used  to  incorporate  [K]bm  similar  to  Chapter  III.  The 
stiffness  coefficients  of  [K]bm  are  evaluated  using  the  analytical  expressions  presented 
in  Chapter  n.  Each  rotating  shaft  is  modeled  as  an  Euler  beam  in  the  lumped  parameter 
model  and  as  a  Timoshenko  beam  in  the  dynamic  finite  element  model  but  the  effect  of 
the  gyroscopic  moment  on  the  shaft  dynamics  is  not  included  The  casing  is  assumed  to 
be  rigid  in  the  lumped  parameter  model  and  flexible  in  the  dynamic  finite  element  model. 
The  rigid  casing  assumption  is  valid  for  many  practical  designs  as  it  may  be  massive  and 
the  rotating  speeds  may  be  sufficiently  low  to  avoid  significant  elastic  deformation.  The 
governing  equations  for  both  discrete  vibration  models  can  be  given  in  the  general 
matrix  form  as 

IM][q(t)}a  +  [C]{q(t)}a  +  [K][q(t)}a=  [f(t)}a  (4.1) 

where  [M],  [C]  and  [K]  are  the  system  mass,  damping  and  stiffness  matrices 
respectively,  and  {q(t)}a  and  {f(t)}a  are  defined  as  the  generalized  alternating 


displacement  and  applied  load  vectors  respectively.  Like  Chapter  III  of  this  report, 
{f}bm  and  bearing  preloads  are  also  assumed  constant  to  ensure  a  time-invariant  [K]bm 
matrix.  Energy  equivalent  viscous  damping  matrix  [C]b  =  a  [K]bm  is  assumed  for  the 
energy  dissipation  mechanism  in  the  bearing  where  o  is  the  Rayleigh  damping 
proportionality  constant.  Other  features  of  the  proposed  theory  given  in  Chapter  III  of 
this  report  are  retained. 

The  specific  objectives  of  this  study  are  to:  (i)  incorporate  the  proposed  bearing 
matrix  [K]bm,  developed  in  Chapter  II,  in  discrete  vibration  models  of  generic  geared 
system  of  Figure  4.1  as  described  by  equation  (4.1)  using  lumped  parameter  and/or 
dynamic  finite  element  method,  (ii)  conduct  overall  system  studies  by  calculating 
eigensolutions  and  forced  harmonic  responses,  (iii)  evaluate  the  effects  of  casing  and 
mount  dynamics  on  the  internal  rotating  system,  and  predict  vibration  transmission 
through  bearings,  (iv)  examine  the  following  3  example  cases:  I.  single-stage  rotor 
system  with  flexible  shaft  supported  by  two  identical  rolling  element  bearings  on  rigid 
casing  and  flexible  mounts  as  shown  in  Figure  4.2,  II.  spur  gear  pair  system  with  motor 
and  load  inertias  attached  to  two  flexible  shafts  supported  by  four  rolling  element 
bearings  on  flexibly  and  rigidly  mounted  rigid  casing  as  shown  in  Figure  4.1,  and  III. 
case  II  with  flexible  casing  and  rigid  mounts,  and  (v)  demonstrate  the  advantages  of  our 
formulation  over  the  existing  vibration  models. 
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4.4  THEORY 

4.4.1  Method  A:  Lumped  Parameter  Model 

4.4.1. 1  Equations  of  Motion 

Both  physical  systems  shown  in  Figures  4.1  and  4.2  are  discretized  using  lumped 
parameter  technique  to  yield  equation  (4.1)  through  the  Lagrange's  equations  of  motion 
[45]. 


a<lw  J 


3Et  3Eu  3Ec 

3qw  3qw 


w  =  1,2,3,... 


(4.2) 


where  ET  and  Ey  is  the  kinetic  and  potential  energies  respectively, 
Ec  =  -j{q]  [C] {q}  is  the  Rayleigh's  dissipation  function  and  Fw  is  the 
generalized  force.  The  total  system  potential  energy  Eu  and  kinetic  energy  Ej  are 
obtained  by  adding  the  energy  of  each  system  component  which  are  derived  in 
subsequent  sections.  Equation  (4.2)  when  applied  to  Figure  4.1  will  lead  to  equation 
(4.1).  It  is  then  rearranged  in  terms  of  the  following  {q}  to  obtain  the  partitioned 
system  matrices 


{q(0}a={{u(t»^{«(,)}^ 


(4.3a) 


[K]  = 


=  DIAG  ■ 

[K  i™)!  ml} 

(4.3b) 

r 

to. 2  to13i 

[K]2, 

TO 22  [K]23 

;  [C]  =  O  [K] 

(4.3c) 

L 

tK]32  [K]33J 
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where  {©(t)}Ra  (subscript  R  indicates  rotor)  consists  of  the  alternating  angular 
displacement  of  all  rotors  including  gears,  motor  and  load;  (m)  and  {1}  are  the 
lumped  mass  and  inertia  row  vectors  respectively;  and  the  operator  DIAG{  )  transform 

the  row  vector  into  a  diagonal  matrix  with  components  of  the  vector  corresponding  to 

T 

the  diagonal  elements.  The  stiffness  sub-matrices  [K]wj  =  fK]^w,  w,j  =  1,  2,  3, 

consist  of  the  appropriate  terms  corresponding  to  the  partitioned  {q(t))a.  For  example, 
T  T  T 

[K]n  couples  { (u(t) }  (u(t) )  )  degrees  of  freedom  to  itself,  and  [K]12  provides 

a  coupling  between  { (u(t)  (u(t)  }T  and  { (BCt))^  (0(t)  }T  degrees  of 

freedom. 

Now  assume  that  the  rotary  inertia  of  each  shaft  lumped  mass  is  negligible.  The 
system  matrices  given  in  equation  (4.3)  arc  therefore  simplified  by  eliminating  {©(t))^ 
degrees  of  freedom  from  {q(t))a  in  equation  (4.3a).  The  resulting  mass  and  stiffness 
matrices  are 


[M]  =  DLAG  {  (m) ^  {m}^  I  (I)TR{I)J} 


(4.4a) 


[K]  = 


[K)11-[K)13[K]-;iK]31  [K]n- 

L  W2rW23[Khl[Khl  TO 22-™ 23^33^32 


(4.4b) 


This  analytical  lumped  parameter  formulation  will  be  used  for  example  cases  I  and  II.  It 
may  be  noted  that  if  the  effect  of  the  rotary  inertias  is  required,  equation  (4.3)  can  be 
used  instead. 


107 


4.4. 1.2  System  Kinetic  Energy 

A  flexible  shaft  of  length  Ls  with  rigid  rotor  of  mass  mR  is  subdivided  into  ns 

number  of  segments  of  equal  length  Le  =  Ls/ns  with  lumped  masses  at  both  ends  of  each 

segment  as  shown  in  Figures  4.3a  and  4.3b  for  a  single-spur  gear  mesh  geared  drive 

and  single-stage  rotor  system  respectively.  Each  lumped  mass  has  3  translational  and  3 

rotational  degrees  of  freedom  as  illustrated  in  Figure  4.3c.  The  total  system  kinetic 

T 

energy  Ej  is  given  by  ^  { 4}  [M]  {  q} .  The  mass  matrix  [M]  is  diagonal  and  consists 
of  the  lumped  masses  and  inertias. 

4.4. 1.3  Shaft  Stiffness  Matrix 

e 

Using  the  direct  stiffness  approach,  the  stiffness  matrix  [K]  s  of  dimension  12 
corresponding  to  the  alternating  displacement  vector  (q(t)  =  (uxj(t),  uxj+1(t),  uyj(t), 

Uyj+i(0,  uzj(t),  uZj+1(t),  0yj(t),  0yj+i(t),  exj(t),  exj+i(t),  e2j(t),  ezj+i(t)}T,  of  a  generic 
shaft  segment  is  given  by 


[K] 


e 

s 


tKu/s  (Ku/s 

fKe/s  tKe/s_ 


(4.5) 


and  the  non-zero  elements  of  [Kuufs,  [  Ku0]®  =  [  K0i|]*T  and  [K^]®  of  dimension 
6  are  given  by  the  following  equations  (4.6),  (4.7)  and  (4.8)  respectively. 


kj  i  =  k?2  =  k33  =  k44  =  -kj2  =  — k2i  =  -*34  =  -k43  =  12  El/L^  ; 
^55  =  ^66  =  “^56  ~  — k^5  =  AE/Le 


(4.6) 
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Bearing 
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Rotor 
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Bearing 

[s+ 
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ns+l 


Flexible  Shaft 


(b) 


x 


A 

y 


►  z 


Figure  4.3  Lumped  parameter  model  of  the  flexible  shafts  with  rigid  rotors.  Each  shaft 
is  supported  by  2  identical  bearings  on  both  ends,  (a)  Internal  rotating 
system  of  a  typical  geared  rotor  system  of  Figure  4.1.  (b)  Internal  rotating 
system  of  a  typical  single-stage  rotor  system  of  Figure  4.2.  (c) 

Discretization  of  a  shaft  segment  and  the  degrees  of  freedom  associated  with 
each  lumped  mass. 
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kl  1  =  -k22  =  k33  =  -k44  =  k12  =  -k2l  =  k34  =  _k43  =  6  El/L  \  (4.7) 

kH=k22=k33  =  k44  =  4El/Le  ; 

k12  =  k21  =  k34  =  k43  =  2  El/L*.  ; 

k55  =  *66  ~  -k56  =  —k65  =  GJ/Le  (4.8) 


where  E  is  the  modulus  of  elasticity,  G  is  the  shear  modulus  of  elasticity,  I  is  the 
moment  of  inertia,  and  J  is  the  polar  moment  of  inertia  of  the  shaft.  The  lumped 
stiffness  matrix  [K]s  corresponding  to  (q(t))sa  =  (uxj(t),  uyj(t),  uzj(t),  9yj(t),  9xj(t), 

T  C 

0zj(t)}  ,  j  =  1,2,3,...  for  the  shaft  is  constructed  by  the  superposition  of  all  [K] 
matrices  and  merging  terms  associated  with  each  degree  of  freedom. 


[K] 


s 


lKu»U 

Kuh  tKeel5  J 


(4.9) 


T 

The  non-zero  elements  of  [Kuu]s,  [K  u0]  s  =  [K  and  [Kg0]s  of  dimension  3(ns+l) 
are  similarly  given  by  the  following  equations  (4.10),  (4.11)  and  (4.12)  respectively. 


kll  ~  k22  ~  k3ns+l,3ns+l  “  k3ns+2,3ns+2  “  ^2  El/Le  ; 

k33  =  k3ns+13ns+l  =  AE/Le  (4.10a) 

kwj+3,wj  =  kwj,W[+3  =  —  ^2  El/Le  ;  wj  =  l,2,4,5,7,8,...,3ns-l  (4.10b) 

kwjj,wjj =  24  El/L  e  ,  Wjj  —  4,5,7,8,10, 1  l,...,3ns— 1  (4.10c) 

kwra+3,wra  —  kwm,wm+3  ~  — AE/Lc  ;  wjjj  =  3,6,9,...,3ns  (4.10d) 

kw,v.wiv  =  2  AE/Le  ;  wjy  =  6,8,12,...,3ns  (4.10e) 
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2 

kll  =  k22  =  -k3ns+13ns+l  =  -k3ns+2,3ns+2  =  6EI/Le  (4.1  la) 

-kwi+3.wI  =  kwI,wI+3  =  6El/Lg  ;  wr  =  1,2, 4,5,7, 8,...,3ns-l  (4.11b) 

kl  1  =  k22  =  k3ns+l,3ns+l  =  k3ns+2^ng+2  =  4  EIA"e  5 

k33  =  k3ns+ 1 ,3ns+ 1  =  GJAe  (4- 1 2a) 

kwI+3lWl  =  kwItwj+3  =  2  El/Le  ;  W!  =  l,2,4,5,7,8,...,3ns-l  (4.12b) 

kWn.Wn  =  8  E^e  ;  w„  =  4,5,7,8,10,11 3ns-l  (4.12c) 

kwm+3,wm  =  kwffl,wm+3  =  “ GJ/Ee  *  WIII  =  3,6,9,.. .,3ns  (4.1 2d) 

kwjy.wjv  =  2  GJ/Le  5  wjv  =  6, 8,12,..., 3ns  (4.12e) 


Note  that  here  [K]s  is  only  for  a  shaft,  of  dimension  6(ns+l).  The  same  formulation 
should  be  applicable  to  a  geared  system  with  multiple  shafts  where  each  shaft  stiffness 
matrix  is  still  given  by  [K]s.  Accordingly,  the  shaft  potential  energy  Eus  is  given  as 

Eus=lt<[K]>{q}!1  <4-13> 


Formulations  for  the  motor  and  load  rotary  inertias,  and  flexible  coupling  torsional 
stiffness  will  be  given  in  the  later  sections  on  example  cases. 

4. 4.1.4  Gear  Mesh  Stiffness  Matrix 

For  the  generic  geared  rotor  system  with  a  spur  gear  pair,  the  driving  and  driven 
shafts  are  coupled  via  a  linear,  time-invariant  spur  gear  mesh  stiffness  kj,  *  kj,(t).  Now 
we  define  the  gear  mesh  coupling  stiffness  matrix  [K]j,  as  follows  where  coupled 


Ill 


torsional  and  translational  motions  of  the  gear  and  pinion  are  defined  by  { uyGa(t), 

UyPa(t),  0zGa(O»  6zPa(0 1 


tK]h  = 


-  k, 


khdG 


khdP 


-  k. 


khdG 


khdP 


khdG 


khdG 


k»4 


khdPdG 


khdp 


khdp 


khdPdG 


khdP 


(4.14) 


where  dc  and  dp  are  the  gear  and  pinion  diameters  respectively;  subscript  G  and  P  refer 
to  gear  and  pinion  respectively.  Here,  the  y  axis  is  taken  to  be  parallel  to  the  gear  mesh 
force  which  is  along  the  pressure  angle  direction  at  the  mesh  point  The  potential  energy 
EUh  in  this  case  is 


E 


UyPa’  0zGa,0zPa}^hj 


yGa 

lyPa 

'zGa 


zPa 


(4.15) 


4.4.1.5  Flexible  Mount  Stiffness  Matrix 

The  flexible  mounts  are  represented  by  a  diagonal  stiffness  matrix  [K]v 

corresponding  to  (q(t))ca  which  consist  of  effective  stiffness  coefficients  kvw,  w  =  x,  y, 
z,  0X,  0y,  0Z.  This  modeling  procedure  assumes  that  no  coupling  exists  between  the 
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casing  rigid  body  degrees  of  freedom  due  to  the  flexibility  of  the  mounts.  Accordingly, 
one  has 


E 


Uv 


1 

2 


{q>^[K]v{q> 


ca 


(4.16) 


4. 4.1. 6  Bearing  Stiffness  Matrix 

The  proposed  bearing  stiffness  matrix  [K]bm  corresponding  to  the  bearing  degrees 
nr 

of  freedom  (5wa(t),Pwa(t)}  ,  w  =  x,  y,  z,  has  been  formulated  in  Chapter  II.  The 
corresponding  potential  energy  EUb  is 


Ub 


waj 

P 

^waj 


(4.17) 


Here,  { Swaj(t),pwaj(t) } T  may  be  expressed  in  terms  of  {qW}^  and  {q(t))ca  through  a 

— > 

coordinate  transformation  for  the  j-th  bearing  located  at  Rj  =  {xj,  yj,  zj)  from  the 
casing  center  of  mass.  For  the  geared  rotor  system,  four  bearings  are  denoted  by  j  =  1 , 
ns+l,  ns+2  and  2ns+2  as  shown  in  Figure  4.3a,  while  for  the  single-stage  rotor  system 
two  bearings  are  denoted  by  j  =  1  and  ns+l  as  shown  in  Figure  4.3b.  Assuming  the 
small  angle  approximation,  {5waj(t),Pwaj(t)}T  for  the  j-th  bearing  is 
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(4.18) 
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4.4.2  Method  B:  Dynamic  Finite  Element  Formulation 

The  dynamic  finite  element  technique  of  incorporating  [K]bm  in  equation  (4.1)  is 
employed  especially  when  the  casing  plate  is  elastically  deformable  over  the  frequency 
range  of  interest,  as  discussed  previously  in  Chapter  III.  In  this  method,  [K]bm  is 
implemented  in  the  dynamic  Fmite  element  model  as  a  generalized  stiffness  matrix.  This 
stiffness  element  provides  the  only  coupling  between  the  flexible  shaft  and  casing  at  the 
appropriate  nodal  points  on  both  components  corresponding  to  the  bearing  locations. 
Additionally,  the  gear  mesh  stiffness  kh  in  Figure  4.1  behaves  like  a  linear  translational 
spring.  Since  the  coupling  between  torsional  and  transverse  shaft  vibrations  is  due  to  kh 
only,  the  corresponding  portion  of  the  stiffness  matrix  is  also  given  by  equation  (4.14). 
Our  finite  element  formulation  uses  conventional  structural  elements  typically  available 
in  commercial  software  programs  [39]  for  the  shaft,  casing  plate  and  mounts  -  this  will 
be  illustrated  for  example  cases  I,  II  and  HI.  Other  features  of  this  method  are  similar  to 
those  discussed  earlier  in  Section  4.4.1. 

4.4.3  Other  Methods 

Alternate  techniques  of  formulating  the  geared  rotor  system  dynamic  problems 
such  as  finite  difference  (which  is  similar  to  method  B),  flexibility,  component  mode 
synthesis  and  transfer  matrix  formulations  are  also  possible.  It  may  be  noted  that 
Berman  [54]  has  philosophically  proposed  application  of  the  component  mode  synthesis 
to  geared  problems.  Other  researchers  [40,41]  have  utilized  the  transfer  matrix  method 
in  the  rotor  dynamic  problems.  However,  such  models  need  to  be  modified  or  extended 
in  order  to  solve  the  overall  dynamic  problems  discussed  in  this  research.  Direct 
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application  of  these  alternate  methods  [40-43,54]  are  beyond  the  scope  of  this 
report  and  are  left  for  further  research. 

4.5  EXAMPLE  CASE  I:  SINGLE-STAGE  ROTOR  SYSTEM  WITH 
RIGID  CASING  AND  FLEXIBLE  MOUNTS 

4.5.1  Vibration  Models 

Consider  the  single-stage  rotor  system  shown  in  Figure  4.2  where  two  axially 
preloaded  identical  deep  groove  ball  bearings  support  a  flexible  rotating  shaft;  see  Table 

4. 1  for  bearing  parameters.  The  shaft  is  assumed  to  be  subjected  to  a  mean  torque  T^,,, 
*  TxsmCO  and  harmonic  excitation  forces  {F^t)  =  F^^QJ  e^1,  Fysa(t)  =  Fysal(Qz) 
ei(Qzt-n/2)  }T  ^ue  tQ  ^  rotating  rotor  mass  unbalance.  The  rotor  is  fixed  to  the  center  of 
the  shaft  which  coincides  with  the  rigid  casing  centroid.  The  shaft  is  assumed  to  be 
transversely  decoupled  from  the  motor  and  load  due  to  the  flexible  torsional  couplings. 

A  lumped  parameter  model  with  ns+2  lumped  masses  is  developed  according  to 

the  theory  given  in  Section  4.4.1.  The  bearing  stiffness  matrix  for  an  axially  preloaded 
ball  bearing  has  significant  stiffness  coefficients  kbxx,  kbyy,  kbzz,  kb0x0x,  kb0y0y,  kbx0y 
and  kby0x  which  depend  on  the  mean  axial  preload.  Since  there  is  no  external  Fj^t) 

force  on  the  system  and  no  coupling  between  the  torsional  and  transverse  motions  of  the 
shaft,  the  dynamics  associated  with  {u2saj(t),uzca(t),02Ra(t))T  degrees  of  freedom  (here 
subscript  j  is  a  dummy  index  to  identify  the  shaft  lumped  masses)  are  decoupled  from 
the  others  and  have  trivial  steady-state  particular  solution.  The  system  matrices  of 
equation  (4)  can  now  be  simplified  by  suppressing  such  degrees  of  freedom.  Also,  it 
can  be  easily  observed  from  equations  (1)  and  (4)  that  2  sets  of  uncoupled  differential 
equations  exist.  One  set  has  {q(t)}a  =  {uxsaj(t),uxca(t),0yRa(t),0yca(t)}T  which  are 
excited  by  Fxsa(t)  and  the  other  set  consists  of  {q(t)}a  =  [uysaj(t),  uyca(t),  0XRa(t), 


115 


Table  4. 1  Design  parameters  for  deep  groove  ball  bearings  used  for  example  cases  ^ 


Load-deflection  exponent  n 

3/2 

Load-deflection  constant  Kn  (N/mn) 

7.5  E9 

Number  of  rolling  element  Z 

12 

Radial  clearance  rL(mm) 

0.005  -  0.02 

Pitch  diameter  (mm) 

38.5 

Aq  (mm)  tt 

0.0625 

Unloaded  contact  angle  0^ 

23° -47° 

t  These  bearings  are  currently  being  used  in  a  NASA  gear  test  facility  [55].  Also  see 
example  cases  II  and  HI. 

tt  Unloaded  distance  between  inner  and  outer  raceway  groove  curvature  centers. 
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0XCa(O}  with  Fysa(t).  The  steady-state  solution  to  these  two  sets  will  always  result  in 
the  synchronous  whirling  motion  of  the  internal  rotating  shaft.  However,  since  the  two 
sets  of  equations  are  similar  and  independent,  only  one  set  is  considered  in  the  following 
analysis.  Hence,  the  problem  reduces  to  a  2  dimensional  vibration  system  with  DOF  = 
ns+4  which  consists  of  a  flexible  shaft  with  rigid  rotor  vibrating  transversely  and 
coupled  to  a  rigid  casing  and  a  flexible  mounts  through  two  bearings.  Accordingly,  [M] 
in  equation  (4.4)  is  rewritten  in  terms  of  the  alternating  displacement  vector  {q(t)}a  = 
{ uxsajW.uxca(t)»0yRa(Oi0yca(t) ) 

[M]  =  DIAG  (m  ^  ,  m  c>  I  yR  » I yc)  (4.19) 

where  m  and  Iy  are  the  lumped  mass  and  inertia  about  the  y  axis  respectively,  and  the 

subscript  e  refers  to  the  shaft  segment.  To  define  [K]  in  equation  (4.4),  the  nonzero 
elements  of  [K]n,  [K]12,  [K]22,  [K]23,  [K]13  and  [K]33  are  given  by  equations  (4.20), 

(4.21),  (4.22),  (4.23),  (4.24)  and  (4.25)  respectively. 


3 

kll  ~  kns+l,ns+l  —  12  El/Le  +  k^xx  ,  kns+2>ns+2  =  kvx  2  kjjXX  ; 

kns+2,l  ~  kns+2,ns+l  —  ^1  jis+2  =  kns+l  ,ns+2  =  — kbxx  (4.20a) 

kwj+l,wj  =  ^wj.wj+l  =  — 12  El/L  c  ;  wj  =  1,2,3, ...,ns  (4.20b) 

=  24  EI/L  e  ;  wjj  =  2,3,...,ns  (4.20c) 

2 

kns/2,l  =-kns/2+2,l  =6EI/Lc  ;  kns+2>2  =  (zj  +  z„s+1)  kbxx  +  2  kbx0y  ; 

k12  =  -Z1  kbxx  “  kbx8y  5  kns+l,2  =  -zns+l  kbxx  -  kbx0y  (4-21) 


k j  j  =  8  El/Le  ;  k22  =  kv6y+2  k^^  (z2  +  z2  s+1)kbxx  +  (zi+zns+1)  kbx0y  (4.22) 


kl,ns/2  =  ^1, ns/2+1  ~  2  El/Le  ; 

k2i  =  -zi  kbxQy  _kbQy9y  ;  k2^,s  =  -Zns+i  kbx0y  _kb0y0y  (4.23) 

2 

kl  1  =  “kns+l,ns  =  6  EI/L  e  +  kbx0y  ;  kns+2j  =  kns+2>ns  =  -kbx0y  , 

kns/2+l, ns/2+1  =  _kns/2+l,ns/2  =  6  El/L \  (4.24a) 

2 

k12  =  ~ kns+l,ns-l  *  -kns/2.ns/2-l  =  kns/2+2jis/2+2  =  6  El/L  e  ;  ns>2  (4.24b) 

2 

kwIU-l ,wUI  =~kwnj-l,win-2  “_kwIn+ns/2,wIu-2+ns/2  =  ^m+nsAwjjj+ns^  -  6  El/L  e  , 

ns>6,  wm  =  3,4,-,ns/2  (4.24c) 

kl  1  =  knStns  =  ^  El/Le  +  kb8y0y  (4.25a) 

kwiV,wiv  =  8  EI/Le  I  ns>2,  W1V  *  2,3,...,ns-l  (4.25b) 

kwy.wy— 1  =  2  E^A^e  ♦  kwy-l,wy  =  2  El/L^  * 

ns>2,  wy  =  2,3,... ,ns  (wy  *  ns/2+l)  (4.25c) 


T 

Due  to  the  symmetric  nature  of  [K],  [K]  .  =  [K] .  .  The  generalized  load  vector 

WJ  JW 

{f(t)}a  =  {0,...,  Fxsa(t),  0,...}  consists  of  only  the  x  component  of  the  rotor  mass 
unbalance  force.  The  other  set  of  linear  governing  equation  may  be  obtained  by 
interchanging  subscripts  x  and  y  in  the  above  formulation.  Our  proposed  vibration 
model  can  also  be  readily  reduced  to  the  conventional  rotor  dynamic  models,  excluding 
the  gyroscopic  moment,  internal  structural  damping,  dissimilar  bearings  and  other 
secondary  effects,  by  retaining  only  kbxx,  kbyy  and  kbzz  in  our  bearing  formulation. 
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4.5.2  Eigensolution 

The  natural  frequencies  C0j  and  modes  <j>j  of  Figure  4.2  with  system  parameters 
given  in  Table  4.2  are  obtained  numerically.  Initially,  we  check  for  coj  convergence  rate 
as  shown  in  Figure  4.4.  through  0)4  and  co5  by  varying  the  number  of  shaft  segments 
ns;  it  may  be  noted  that  <»i,  0)2  and  CO3  behave  like  0)4,  and  tog  trend  is  very  similar  to 
(05.  Beyond  ns  =  4,  there  are  no  noticeable  changes  in  the  natural  frequencies.  These 
predictions  are  now  compared  with  two  simple  lumped  parameter  vibration  models  and  a 
dynamic  finite  element  model  (FEM)  in  Table  4.3.  One  of  the  simple  models  excludes 
the  effect  of  casing  and  flexible  mounts,  and  both  employ  a  conventional  bearing  model 
with  only  kbxx  incorporated  as  described  earlier  in  Section  4.5.1.  The  dynamic  finite 
element  formulation  models  the  flexible  casing  constructed  with  four  noded  quadrilateral 
plate  elements  -  a  description  of  the  methodology  will  be  presented  later  in  example  case 
III.  Our  proposed  theory  and  FEM  are  in  excellent  agreement  with  each  other  as  the 
deviations  are  within  ±4%  for  first  6  natural  frequencies.  The  corresponding  mode 
shapes  <J>j  are  shown  in  Figure  4.5.  Note  that  for  each  mode  in  the  x-z  plane  given  here, 

there  is  a  complementary  mode  in  the  y-z  plane  which  will  have  the  same  natural 
frequency  if  kvy=kvx,  kv0j(=kv0y  and  Ixc=IyC.  Using  the  exactly  same  parameters, 

results  obtained  from  two  simple  models  are  also  compared  in  Table  4.3.  The  simple 
model  with  casing  and  mount  dynamics  included  predicts  lower  natural  frequencies  as 
compared  to  our  prediction  except  for  W5;  here  the  symbol  a  implies  estimation  using  a 
simple  formulation.  This  model  gives  a  rough  approximation  for  most  of  the  modes,  as 
evident  from  Table  4.3.  Further,  the  simple  model  without  casing  and  mounts  can  not 
predict  all  the  modes  below  2  kHz.  In  addition,  to 3  shifts  while  other  modes  remain 
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Table  4.2  System  parameters  of  example  case  I:  single-stage  rotor  system  shown  in 
Figure  4.2 


Rotor  mass  hir  (kg)  and  inertia  IyR  (kgm2) 

5.0,  1.24E-2 

Shaft  mass  mg  (kg)  and  length  Ls  (m) 

1.03,  4.2E-1 

Casing  mass  trie  (kg)  and  inertia  Iyc  (kgm2) 

148,  4.27 

Shaft  flexural  rigidity  El  (Nm2) 

1.626E3 

Bearing  axial  preload  Fzbm  (N) 

238 

^bxx’  ^bx9y>  ^b0y0y  ^ 

1.37E8,  -1.96E6,  4.27E4 

^VX>  ^V0y 

4.0E8,  1.8E7 

t  Computed  kbwj  as  proposed  in  Chapter  II.  Only  relevant  bearing  coefficients  are 
listed. 


Table  4.3  Comparison  of  natural  frequencies  in  Hz  for  example  case  I 


Proposed 

Theory 

Simple  Models 

Mode 

FEM 

A,% 

without  casing 

with  casing 

a2% 

1 

97.2 

97.2 

0.0 

69.5 

69.4 

28.6 

2 

262 

262 

0.0 

not  predicted 

262 

0.0 

3 

309 

309 

0.0 

297 

294 

4.9 

4 

332 

334 

0.6 

not  predicted 

329 

1.5 

5 

1405 

1365 

-2.9 

1442 

1443 

-5.7 

6 

1755 

1690 

-3.8 

1452 

1453 

14.0 

A[%  =  100  x  (FEM  -  Proposed  Theory)  / FEM 
A2%  =  100  X  (FEM  -  Simple  Model)  /FEM 


3 


Number  of  Shaft  Segments,  n  s 


Figure  4.4  Natural  frequencies  0)4  and  0)5  versus  the  number  of  shaft  segments, 
example  case  I. 
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Figure  4.5  First  6  mode  shapes  of  example  case  I:  single-stage  rotor  system  including 
rigid  casing  and  flexible  mounts. 
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nearly  the  same.  These  observations  are  consistent  with  the  findings  reported  earlier 
[51-53]. 

We  examine  the  effect  of  casing  mass  ir^  (and  inertia  Iyc  «  rr^)  while  retaining 
other  parameters  of  Table  4.2.  The  fundamental  frequency  (Oj  (first  transverse  mode  of 
the  shaft-rotor  system)  and  <03  (second  transverse  shaft-rotor  mode)  are  quite  insensitive 
to  mj..  On  the  other  hand,  CO5  and  cog  are  proportional  to  (mc)'^2  for  a  light-weight 
casing  given  by  mc/(ms+mR)  <  1,  but  are  invariant  for  mc/(ms+mR)>  1  as  shown  in 
Figure  4.6a.  The  converse  is  seen  for  £1)2  and  004  which  are  nearly  constant  for 
mc/(ms+mR)  <  1,  but  are  proportional  to  (n^)'1/2  for  1^/(1115+0^)  >  1.  These  results 
imply  that  the  single  degree  of  freedom  approximation  can  be  applied  to  CO5  and  0)6  in 

the  region  where  mc/(ms+mR)<  1,  and  similarly  to  0)2  and  oa4  in  the  other  region. 

2 

Next,  we  vary  mount  stiffness  coefficients  kvx  (and  *  L  s  kvx)  to  simulate 
the  effect  of  mount  flexibility  on  the  system.  Here  <05  and  tag  are  found  to  be  unaffected 
by  kvx  due  to  the  fact  that  the  corresponding  modes  are  predominantly  shaft  bending 
motion  type  as  shown  in  Figure  4.5.  Figure  4.6b  indicates  that  oil  and  01 3  are 
proportional  to  (k^)1^2  for  soft  mounts  (kvx/kbxx  «  1)  but  are  insensitive  to  kvx  for 
stiff  mounts  (kvx/kbxx  »  1).  Conversely,  0)2  and  0)4  are  proportional  to  (kvx)'^2  for 
stiff  mounts  but  remain  nearly  invariant  for  soft  mounts.  Similarly,  each  mode  in  its 
appropriate  region  can  be  assumed  to  behave  like  a  single  degree  of  freedom  system.  In 
Figure  4.6b,  unlike  Figure  4.6a  where  only  one  transition  point  at  m^ms+mR  exists,  2 
transition  points  are  found.  These  are  located  below  and  above  kvx/kbxx  =  1,  and 
increase  with  increasing  shaft  and  bearing  stiffnesses. 

Now,  we  examine  the  range  of  bearing  preload  Fzbm  from  100  to  5000N  on  coj.  It 
is  seen  that  coj  increases  with  increasing  Fzbm  due  to  an  increase  in  the  magnitudes  of  the 


Normalized  Mount  Stiffness,  kvx  /kbxx 
(b) 


Figure  4.6  Effect  of  casing  mass  normalized  with  respect  to  shaft  and  rotor  masses, 
and  mount  stiffness  kvx  normalized  with  respect  to  kj,xx  on  system  natural 
frequencies  0);  for  example  case  I. 
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bearing  stiffness  coefficients  k^wj  as  shown  in  Table  4.4  for  the  two  extreme  Fzbm 
values;  this  observation  is  similar  to  those  found  earlier  in  the  bearing  system  studies  of 
Chapter  III.  Note  that  o>2  *s  the  most  insensitive  to  Fzbm,  and  0)3  and  0)4  show  a 
moderate  variation.  But  g>i  increases  by  about  45%  due  to  larger  bearing  motions  as 
compared  to  0)2,  0)3  and  CO4  modes.  The  remaining  two  modes  of  interest  are  also 
affected  significantly  as  large  motions  across  the  bearings  are  again  found  in  these 
modes.  Unlike  the  trends  associated  with  variations  of  n^.  and  kvx,  here  each  coj 
approaches  an  upper  bound  solution,  as  shown  in  Figure  4.7  for  toj,  if  Fzbm  is 
sufficiently  high  to  stabilize  magnitudes  of  kbwj  terms. 


Table  4.4  Effect  of  bearing  preload  Fzbm  on  tOj  (Hz)  for  example  case  I 


Preload  Fzbm  (N)  ©!  C02  w3  ©4  o)6 


100 

89.4 

262 

305 

331 

1366 

1706 

5000 

130 

263 

321 

346 

1855 

1968 
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4.5.3  Transmissibility 

The  forced  harmonic  response  due  to  the  rotating  rotor  mass  unbalance  force 
Fxsa(t)  =  Fxsai(Dz)  e'^zt  is  obtained  using  the  dynamic  stiffness  technique  given  earlier 
in  Chapter  III.  Also,  recall  from  Chapter  III  the  definition  of  sinusoidal  load 
transmissibility  functions  R(QZ)  between  two  arbitrary  locations  I  and  II  with  load 
vectors  fwia(t)  and  fjIla(t)  as  R(QZ)  =  |fwIa(Qz)  l/lfin  a(£lz)|  where  in  this  case 

fjllaW  =  Fxsa(t).  The  bearing  force  RFxba,Fxsa<^z)  and  moment  RMyba,Fxsa(nz) 

transmissibility  spectra  are  compared  in  Figure  4.8  for  our  formulation  and  simple 
model,  given  the  same  system  of  Table  4.2,  both  models  differ  only  in  [K]bm 

formulation.  Note  that  the  results  presented  are  for  only  one  bearing  due  to  symmetry  of 
the  problem.  We  observe  that  the  simple  model  predicts  lower  RFxba,Fxsa(^z)  than 

proposed  model  except  in  the  vicinity  of  a>j.  On  the  other  hand,  RMyba,Fxsa(^z)  *s  not 

predicted  at  all  by  the  simple  model.  This  component  of  the  bearing  transmissibility 
spectra  is  primarily  due  to  coefficients  kbx0y  and  kb0y0y  which  are  obviously  not 

included  in  the  simple  model.  Our  predictions  also  show  that  higher  modes  over  250- 
350  Hz  contribute  significantly  to  the  bearing  transmissibility  which  are  not  seen  in  the 
force  transmissibility  spectra  yielded  by  the  simple  model. 

Now  we  compare  the  mount  load  transmissibility  in  Figure  4.9.  We  observe  that 
the  trends  of  mount  force  transmissibility  RFxva,Fxsa(^z)  spectra  predicted  by  both 

models  are  similar,  although  the  simple  model  predicts  slightly  lower  amplitudes  above 
©i.  The  simple  model  again  can  not  predict  the  moment  transmission  RMyva,Fxsa(^z) 

across  the  mounts.  Unlike  the  bearing  transmissibility  spectra,  here  the  effects  of  $3 
and  4>4  on  p  (Qz)  am  more  significant  than  <|>i  and  <j>2-  This  is  mainly  due  to 

yva*  xsa 

the  large  angular  motions  on  the  casing  for  <£3  and  $4. 


Figure  4.8  Bearing  transmissibility  for  example  case  I.  (a)  Force  transmissibility 
Re  p  («z)-  (b)  Moment  transmissibility  Rm  k  ,f  (&z)  which  is  not 

xbi*  xm  xsa 

predicted  by  simple  model. 

JS 


Moment  Transmissibility  (m)  Force  Transmissibility 
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(a) 


(b) 


Figure  4.9  Mount  transmissibility  for  example  case  I.  (a)  Force  transmissibility 
RFxva.Fx»a(fiz)-  (b)  Moment  transmissibility  RMyva,FXSIW- 


Mount  Transmissibility  Bearing  Transmissibility 
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(a) 


(b) 

Figure  4.10  Effect  of  mount  stiffness  kvx  on  the  force  transmissibilities  for  example 
case  I.  (a)  Bearing,  (b)  Mount. 


Mount  Transmissibility  Bearing  Transmissibility 
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(a) 


Figure  4.11  Effect  of  casing  mass  m^.  on  the  force  transmissibilities  for  example  case  I. 
(a)  Bearing,  (b)  Mount. 
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Next,  we  investigate  the  effects  of  and  kvx  on  the  bearing  and  mount 
transmissibilities.  Figure  4.10a  shows  that  essentially  the  same  bearing  transmissibilities 
are  obtained  for  three  kvx  values.  In  the  case  of  the  mount  transmissibility  as  shown  in 
Figure  4.10b,  the  amplitudes  in  general  decrease  with  lower  ky*  except  near  resonances. 
The  casing  mass  also  does  not  affect  the  bearing  transmissibility  significantly  as 
shown  in  Figure  4.1  la.  The  mount  transmissibility  amplitudes  increase  with  lower 
as  shown  in  Figure  4.11b  except  near  resonances  again.  It  may  be  noted  that  similar 
trends  are  observed  in  the  bearing  and  mount  moment  transmissibility  spectra. 

4.6  EXAMPLE  CASE  II:  GEARED  ROTOR  SYSTEM  WITH  RIGID 
CASING  AND  FLEXIBLE  MOUNTS 

4.6.1  Bearing  Analysis 

Now  consider  the  geared  rotor  system  in  Figure  4.1  with  flexibly  mounted  rigid 
casing  whose  bearing  and  other  system  parameters  are  given  in  Tables  4.1  and  4.5 
respectively.  Four  ball  bearings  which  support  the  two  shafts  are  subjected  to  identical 
mean  axial  displacement  S^.  The  spur  gear  pair  drive  is  driven  by  a  mean  torque  T^ 
*  Tzsm(t)  which  also  generates  mean  radial  bearing  force  Fyi,mj  and  moment  Mxt,mj. 
The  stiffness  matrix  [K]bm  for  each  bearing  under  these  loads  has  significant 
coefficients  kjjXX,  k^yy,  kj,zz,  kb0x0x»  kjj0y0y,  RbxOy*  ^by0x>  ^byz  ^bz8x*  A  set  of 
governing  nonlinear  algebraic  equations,  consisting  of  3  bearing  load-displacement 
relations  for  each  bearing  as  given  in  Chapter  II  and  from  the  shaft  bending  theory  is 
given  by 


^xbmwj  <k})  2EI(2PxmWj+PxmWjj+3(8ymWj— 5ymWjj)/Ls)/Ls  ; 

Mxbmwn=-MxbinWl  ;  for  w„=ns+l  if  wj-1,  and  wn=2ns+2  if  wj=ns+2. 


(4.26) 
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Table  4.5 


System  parameters  of  example  case  II:  geared  rotor  system  shown  in  Figure 
4.1  + 


Gear  and  pinion  masses  mP  =  rriQ  (kg) 

0.5 

Gear  and  pinion  rotary  inertias  IyG=Iyp.IzG=IzP  (k gm2) 

1.5E-4,  3.0E-4 

Number  of  gear/pinion  teeth 

28 

Gear/Pinion  Pressure  Angle  (degrees) 

20° 

Shaft  mass  ms  (kg)  and  length  Ls  (m) 

2.8,  2.54E-1 

Casing  mass  m*.  (kg)  rotary  inertias  1^,1^  (kgm2) 

77,  1.5,  1.9 

Shaft  flexural  rigidity  El  (Nm2) 

1.25E4 

Effective  Torsional  Stiffness  kj^kj,  (Nm/rad) 

6.05E3,  2.45E4 

Motor  and  Load  inertia  1^ 

1.00E-1,  3.35 

Gear  mesh  stiffness  kf,  (N/m) 

1.0E8 

Mean  axial  bearing  displacement  (N) 

6.00E-4 

Mean  input  torque  Tzbm  (Nm) 

72 

kbyy’  kby8*>  kb0x0x  ^ 

9.7E8,  6.0E5,  1.0E6 

kVy,  kvQx,  kv0z 

9.0E8,  1.4E7,  2.4E7 

t  The  internal  rotating  spur  gear  pair  is  currently  being  used  in  a  NASA  gear  test 
facility  [55]. 

ft  Only  relevant  bearing  stiffness  coefficients  are  tabulated. 
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Due  to  the  physical  symmetry,  we  assume  8ymWl  =  8ymWn  and  Pxmwj  =  -Pxmwn  which 
simplify  the  nonlinear  algebraic  equation  set  for  each  bearing  to 

Mxbm  =  Ls  /  (8  <3q)  -  2  El  Pxm  /  Ls  ;  Fybm  =  T^m/d^  (4.27) 

where  Mxbm  and  Fbym  are  expressed  in  terms  of  a  constant  5zm>  and  variables  8ym  and 
Pxm.  Solution  to  these  2  equations  may  then  be  used  to  compute  kbij(8zm,8yrn,Pxm) 
directly. 

4.6.2  Vibration  Models 

A  lumped  parameter  model  of  Figure  4.1  is  developed  in  accordance  with  Section 
4.4.1.  Each  shaft  is  divided  into  ns  number  of  segments  as  shown  in  Figure  4.3a.  The 
system  matrices  of  equation  (4.4)  are  simplified  by  neglecting  the  longitudinal  shaft 
motion  and  casing  rigid  body  degrees  of  freedom  uzca(t),  uxca(t)  and  0yca(t).  Both 
driving  and  driven  shafts  are  assumed  to  be  coupled  to  the  motor  and  load  respectively 
through  flexible  torsional  couplings.  Therefore,  only  the  motor  and  load  1^  rotary 
inertias  are  considered,  and  the  exterior  portion  of  the  shafts  beyond  the  flexible 
couplings  are  modeled  as  purely  torsional  stiffness  elements  which  are  then  combined 
with  the  flexible  coupling  stiffnesses.  The  system  is  excited  by  the  static  transmission 
error  e(t)  at  the  mesh  point  which  generates  gear  and  pinion  force  Fb(t)  =  kh  e(t),  parallel 
to  the  line  of  action,  and  torque  Tj,g(0  =  l/2d<jkhe(t)  or  T|,p(t)  =  l/2dpkbe(t),  about  the 
axial  z-axis.  The  mass  matrix  [M]  in  equation  (4.4a)  in  terms  of  (q(t))  =  {uysaj(t), 
uycaW>  ®xGa(0>  ®zGaW’  ®xPaW*  ®zPaW>  ®zMaW*  ®zLaW’  ®xcaW*  ®icaW^ 
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[M]  —  DIAG  {mscj>mxc.IxG>IzG<IxP>IzP>IzM>IzL»Ixc»^zc  }  (4.28) 

The  matrices  [K]u,  [K]l2,  [K]13,  [K]22,  [K]23  and  [K]33  define  the  symmetric  [K]  of 
dimension  2ns+  1 1  in  equation  (4.4b).  Nonzero  elements  of  [K]n  of  dimension  2ns+3 

are 


kns+2,ns+2  ~  k2ns+2,2ns+2  -  kll  “  kns+l,ns+l  ”  12  El/L  e  +  kbyy  (4.29a) 

kwj+ 1 , wj  =  W{+  i  =  ~  1 2  El/  L  e  . 

*cwI+ns+2,wj+ns+l  =  kWj+ns+l,wI+ns+2  =  El/L  e  ;  Wj  =  1,2,3,. ..,2ns  (4.29b) 

3  3 

kwn,wn  =  24  El/L  e  ;  k  wn+ns+  j iWn+ns+ 1  =  24  El/L  e  , 

wji  =  2,3,...,ns/2,ns/2+2,...,ns  (4.29c) 

3  3 

kns/2+l^s/2+ls=24El/Le  +  kh  ;  k3n^+2,3ns/2+2  =  24  EV L  e  +  kh  » 

kns/2+l,3ns/2+2  =  k3ns/2+2,ns/2+l  * -kh  '*  k2ns+3,2ns+3  =  4kbyy  +  kvy  ;  (4.29d) 

kwffl,2ns4  3  -  k2ns+3,wm  =  -kbyy  *  wm  =  1,  ns+l,  ns+2,  2ns+2  (4.29e) 


The  nonzero  elements  of  [K]j2  of  2ns+3  rows  and  8  columns  are 


£ 

k3ns/2-l,l  -  ~k3ns/2+3,l  =  kns/2,3  =  -kns/2+2,3  =  6  EI/L  e  '*  k2ns+3,7  =  4kby0x  i 
“kns/2+l,2  =  k3ns/2+2^  =  kh/(2dG)  i  kns/2+l,4  =  ~k3ns/2+2,4  =  khA2dp)  (4-30a) 
kwm-7=-kby0x+Zwm  kbyy  »  ^wui^=-Xwiii  kbyy  »  wiil=  ^  ’ns+^ ,ns+2,2ns+2  (4.30b) 


The  nonzero  elements  of  [K]  13  of  2ns+3  rows  and  2ns  columns  are 
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2 

k  1 1  -  kns+2^s+l  -  6  El/L  e  +  kbyQx  ;  k„s+1>ns  =  k2ns+2,2ns  =  El/ Lg  +  kbyQx  ; 

k3ns/2+2,3ns/2+l  =  _k3 115/2+1,31)5/2  =  kns/2+l, ns/2+1  ~  _kns/2+ 1*5/2 =  EI/L  e  ; 
k2ns+3,l  =  k2ns+3,ns  =  k2ns+3,ns+l  =  k2ns+3,2ns  =  “kby9x  (4.31a) 

k12  =  — kns+l,ns-l  =  kns/2, ns/2-1  =  kns/2+2,i>s/2+2  =  kns+2^s+2  =  6  El/L^  ; 

~k2ny+2,2nsrl  =  “k3ns/2+l,3ns/2-l  =  k3ns/2+3,3ns/2+2  ~  6  El/L  c  ;  ns>2  (4.3 lb) 

"^wjv.wj y- 1  =kwj v,wiv+ 1  =kWlv+ns/2+ 1 , wIvns/2+ 1  =_kwIV+ns+l  .wjy-  l+ns^  El/  L  2  ; 

K  wiv+ns+ 1  ,wiv+ 1 +ns~~k wIV+ns/2+ 1  ,wIV+ns/2-  l=kwrv+3ns/2+2,wrv+3ns/2+ 1 =6  El/  L  e  ; 

2 

"kwiv+3ns/2+2,wiv+3ns/2-l  ~  ^  El/Le  ;  nsS6  ;  wjy  =  2,3,...,ns/2-l  (4.31c) 

where  nonzero  elements  given  by  equations  (4.31b)  and  (4.31c)  exist  only  if  ns>2  and 
ns>6  respectively.  The  nonzero  elements  of  symmetric  [Kfo  of  dimension  8  are 

ku=k33  =  8El/Le  ;  k22  =  kTM  +  khd2/4  ;  k25  = -kTM  ;  k^  =  -k^  ; 

k-44  -  kTL+  khdp/4  ;  k24  =  1<42  = -Mg<W4  ;  k55  =  kTM  ;  kft-fcn.  ; 
k77  =  4kb0x0x  +  £  (zj)2  kbyy  +  kv0x  ’»  k88  =  kv9z  +  E  (Xj)2  kbyy  (4.32) 

where  k-pL  and  kj^  are  the  effective  torsional  stiffnesses  due  to  shaft  and  torsional 
coupling  respectively  at  load  and  motor  ends,  and  the  summation  2  is  over  the  4  ball 
bearings  denoted  by  j  =  1,  ns+l,  ns+2,  2ns+2.  The  nonzero  elements  of  [K]23  of  8 
rows  and  2ns  columns  are 

*1*5/2  =  kl,  ns/2+1  =  k33ns/2  =  k3,3n<j/2+l  =  2  El/Le 
k7,wv  =  +  Zwjjj  kby9x  ~  kb9x9x  !  k8,wv  =  “xwm  kby9x  J 


(4.33a) 
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wm=l(wv=l),ns+l(wv=ns),ns+2(wv=ns+l),2ns+2(wv=2ns)  (4.33b) 

Finally,  the  nonzero  elements  of  symmetric  [K^  of  dimension  2ns  are 

^1 1  =  ^ns,ns  =  ^ns+Ms+l  =  ^2ns,2ns  =  ^  El/Le  +  ^b6y9y  (4.34a) 

vi»wvi  =  ^ns+wvi'ns+wvi  =  ^  ^A"e  »  "s>2.  WVI  =  2,3,.. .,ns  (4.34b) 

^wyn.wvir1  =  ^wvir1*wvn  =  kns+wvn-ns+wvrr1  =  2  El/Le  ; 

^ns+wyn-l iis+wyjj =  2  El/Le  ,  ns>2,  wyjj  =  2,3,...,ns  (wyjj  ^  ns/2+l)  (4.34c) 

As  discussed  earlier,  our  proposed  vibration  model  again  differs  from  the  conventional 
spur  gear  pair  dynamic  models  essentially  due  to  [K]bm  formulation.  The  simple  model 
with  casing  and  mount  dynamics  can  be  obtained  from  our  proposed  model  by  retaining 
only  kjjXX,  kbyy  and  kj,zz  in  [K]^. 


4.6.3  Eigensolution 

The  natural  frequencies  coj  and  modes  <j>j  of  Figure  4.1  with  system  parameters 
given  in  Tables  4.1  and  4.5  are  computed  using  the  proposed  theory,  a  dynamic  finite 
element  model,  and  the  simple  theory  with  and  without  casing  and  mounts.  The  FEM 
model  includes  casing  flexibility  -  this  will  be  described  later  in  example  case  II.  Our 
proposed  theory  differs  from  FEM  by  less  than  ±4%  for  the  first  six  ©j  and  is  within 
+10%  for  ©7  and  ©g  as  shown  in  Table  4.6a.  The  corresponding  natural  modes  are 
described  in  Table  4.6b.  Simple  models  in  general  predict  lower  natural  frequencies  and 
deviate  substantially  except  for  ©j  and  ©3.  The  simple  model  without  any  casing 
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Table  4.6  Eigensolution  of  example  case  II 


a.  Comparison  of  natural  frequencies  C0j  in  Hz 


Proposed 
Mode  <{>j  Theory 

FEM 

A,% 

Simple  Models 

without  casing  with  casing 

a2% 

1 

0 

0 

0.0 

0 

0 

0.0 

2 

31.0 

31.0 

0.0 

29.8 

29.1 

6.1 

3 

480 

480 

0.0 

not  predicted 

480 

0.0 

4 

526 

525 

-0.2 

not  predicted 

510 

2.9 

5 

563 

564 

0.2 

not  predicted 

560 

0.7 

6 

790 

760 

-3.9 

666 

650 

14.5 

7 

1147 

1093 

-4.9 

837 

809 

26.0 

8 

1208 

1100 

-9.8 

694 

681 

38.1 

Ai%  =  100  x  (FEM 
A2%  =  100  x  (FEM 

-  Proposed  Theory)  /  FEM 

-  Simple  Model)  /  FEM 

b.  Summary  of  mode  shapes 

Mode  <(>j 

Description 

1 

Motor- gear-pinion- load  system  rigid  body  torsional  motion 

2 

Motor-gear-pinion-load  system  torsional  motion 

3 

Casing  rigid  body  rotational  motion  0xca 

4 

Casing  rigid  body  translational  motion  uyca 

5 

Casing  rigid  body  torsional  motion  0^ 

6 

First  shaft  coupled  transverse-torsional  motion 

7 

Second  shaft  coupled  transverse-torsional  motion 

8 

Shaft  transverse  motion 

dynamics  as  outlined  here  is  in  excellent  agreement  with  the  results  reported  by 
Kahraman  et  al.  [15].  This  provides  further  validation  for  our  theory. 

The  effect  of  casing  mass  ir^  (and  inertias  Ixc,  1^  n^.)  is  shown  in  Figure  4.12a 

for  0)3,  6)5  and  0)7.  Here  ©2  and  ©g  are  unaffected  by  tr^.  On  the  other  hand,  ©3  is 
always  proportional  to  ( 2  except  when  m^^s+mG+mp)  «  1.  This 
proportionality  feature,  similar  to  the  single  degree  of  freedom  theory,  is  also  seen  for 
©7  and  ©g  for  a  lightweight  casing,  and  ©4  and  ©5  for  a  massive  casing.  It  may  be 
noted  that  the  trends  of  ©g  and  ©4  are  similar  to  ©7  and  ©5  respectively  although  they 
are  not  included  in  Figure  4. 12a.  Unlike  example  case  I  where  the  transition  point  is  at 
m^ms+mG+mp)  -  1,  two  transition  points  are  observed  here  due  to  the  effects  of  the 
gear,  pinion  and  casing  rotary  inertias. 

The  effect  of  mount  stiffness  kvy  (and  kVQx,  kvgz  °e  kvy)  on  ©j  is  given  in  Figure 
4. 12b.  We  again  find  that  ©2  and  ©$  are  not  affected  by  a  variation  in  kvy.  Also,  two 
transition  points  are  found  which  separate  a  region  of  constant  ©j  from  the  region  where 
©j «  (kvy)1/2  for  some  of  the  natural  frequencies.  Here  ©g  is  similar  to  ©7  which  is 
proportional  to  (kyy)1/2  for  stiff  mounts  but  is  insensitive  to  kyy  for  softer  mounts.  The 
converse  is  seen  for  ©2  with  a  lower  transition  point  On  the  other  hand,  ©5  «  (kvy)1/2 
in  the  region  between  the  two  transition  points,  and  is  constant  elsewhere.  Finally,  ©3 
is  seen  to  be  always  increasing  with  kyy. 

The  effect  of  mean  bearing  displacement  8^  is  summarized  in  Table  4.7  for  two 
extreme  $zm  values.  Only  the  natural  frequencies  associated  with  the  shaft  transverse 
and/or  torsional  motion  are  sensitive  to  8zm  or  preloads  Fzbm(8zm  ).  The  rate  of 
increase  for  ©j  is  high  for  low  8^,  but  the  rate  decreases  with  a  higher  8zm  as  observed 
earlier  in  example  case  I  due  to  the  stabilization  of  kyj. 


Natural  Frequency  (kHz)  Natural  Frequency  (kHz) 


139 


(a) 


(b) 

Figure  4.12  Effect  of  casing  mass  (normalized  with  respect  to  shaft  and  gear 
masses)  and  mount  stiffness  kvx  (normalized  with  respect  to  kj>xx)  on  the 
system  natural  frequencies  a)j  for  example  case  II. 
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Table  4.7  Example  case  II:  effect  of  mean  bearing  displacement  on  coj  (Hz) 


5/bm  (mm) 

CO! 

C02 

6)3 

(04 

<*5 

<»6 

co7 

o>8 

29.9 

520 

562 

764 

852 

852 

1.00 

480 

526 

563 

791 

1180 

1215 

4.6.4  Transmissibility  Spectra 

The  response  due  to  the  static  transmission  error  excitation  e(t)  at  the  mesh  point  is 
computed  using  the  same  dynamic  stiffness  technique  used  for  example  case  I.  Here, 
only  the  fundamental  harmonic  of  e(t)  at  mesh  frequency  Wj,  is  included.  All 

uansmissibility  functions  are  normalized  with  respect  to  the  magnitude  of  the  gear  mesh 
force  Fh(t)  =  kh  e(t).  Figure  4.13a  and  4.13b  compare  the  bearing  force  RFyba,Fh(wh) 

and  moment  transmissibility  spectra  as  predicted  by  our  formulation  and 

simple  model.  It  may  be  noted  that  the  transmissibility  spectra  for  all  four  bearings  are 

similar.  Although  simple  model  compares  reasonably  well  with  our  proposed  theory  for 
RFyba,Fh<“h).  R  is  not  capable  of  predicting  ^Mxba,Fh(£0h)>  as  also  seen  previously.  In 

the  case  of  the  mount  transmissibility  shown  in  Figure  4.13c,  only  the  net  moment 
Mzva(t)  is  transmitted.  The  resultant  vertical  force  Fyva(t)  and  moment  Mxva(t)  are 
negligible  due  to  the  force  and  moment  cancellations  at  the  mount  feet.  The  Fourier 
spectrum  of  the  normalized  dynamic  transmission  error  p(t)  =  [yP(t)~ 
yG(t)+(^p0p(t)/2)-(dQ0G(t)/2)]/e(t)  is  shown  in  Figure  4.13d.  Only  two  modes,  <>2 
which  is  predominantly  torsional  vibration  of  the  shafts  and  <j>g  which  is  a  coupled 


Mount  Transmissibility  (m)  Force  Transmissibility 
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Figure  4.13  Comparison  of  transmissibility  and  normalized  dynamic  transmission  error 
spectra  for  example  case  II.  (a)  Bearing  force  transmissibility 
RFyb„Fh(°)h)-  (b)  Bearing  moment  transmissibility  RMxba,Fh(Wh)  which 
is  not  predicted  by  simple  model,  (c)  Mount  transmissibility 
RMzb»*Fh^<0h)-  ^  Normalized  dynamic  transmission  error  pfcoj,). 
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Figure  4.14  Mount  transmissibility  studies  for  example  case  II.  Effect  of  (a)  casing 
mass  mg,  (b)  mount  stiffness  kvx  and  (c)  mean  bearing  displacement  8zm 

on  RMzb„Fh(«h)- 


Moment  Transmissibility  (m)  Force  Transmissibility 


Figure  4.15  Effect  of  mean  bearing  displacement  8zm  on  (a)  bearing  force 
transmissibility  RFybt,Fh(a)b)»  (b)  bearing  moment  transmissibility 

^Mxb«>Fh^0>h^  ^  normalized  dynamic  transmission  error  p(coh)  for 
example  case  II. 


r  ■« 

* 
'  V 
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transverse-torsional  mode,  affect  p(coj,)  significantly.  The  sixth  mode  produces  a  higher 
p(cc>h)  amplitude  at  G>h  than  <t>2-  Additionally,  the  simple  theory  underestimates  p((Oj,)  as 
compared  to  the  proposed  model. 

The  effects  of  m*.,  kvx  and  8^  on  the  transmissibility  spectra  and  p(o)j,)  are  given 
in  Figures  4.14  and  4.15.  Mount  transmissibility  Rm-^Fj/0^)  is  lower  for  a  heavier 

casing,  more  compliant  mount  and  higher  8zm  as  shown  in  Figure  4.14.  In  addition, 
the  effects  of  nr^  and  kvx  on  RMzva,Fh(coh)  arc  more  prominent  than  those  of  8zm.  The 

bearing  transmissibility  spectra  and  p(coh)  are  affected  only  by  8zm  as  shown  in  Figure 
4. 15,  and  are  virtually  insensitive  to  and  kvx  for  the  parameters  given  in  Tables  4.1 
and  4.5.  Vibratory  force  and  moment  transmission  through  the  bearing  are  slightly 
lower  for  a  larger  8zm  due  to  the  additional  constraint  provided  by  higher  kbij  values  on 
the  rotating  shafts.  Conversely,  normalized  dynamic  transmission  error  p(iOj,)  increases 
if  a  larger  bearing  preload  is  specified. 

4.7  EXAMPLE  CASE  III:  GEARED  ROTOR  SYSTEM  WITH  RIGIDLY 
MOUNTED  FLEXIBLE  CASING 

4.7.1  Physical  Setup 

The  final  example  case  examines  the  NASA  gear  test  facility  as  shown  in  Figure 
4. 16  [55].  The  system  parameters  are  equivalent  to  those  in  Tables  4.1  and  4.5  except 
for  the  flexible  steel  casing  of  approximate  dimensions  0.33  x  0.28  x  0.25  m  and  plate 
thickness  of 0.006m.  Rigid  mount  feet  attach  the  four  comers  of  the  bottom  casing  plate 
to  a  massive  foundation.  High  precision  gear  and  pinion  are  used  which  are  identical 
with  0.006m  facewidth,  0.089m  diameter  and  1:1  ratio.  Four  axially  preloaded  high 
precision  deep  groove  ball  bearings  are  being  used  to  support  0.03m  diameter  shafts  of 
length  0.254m  on  the  flexible  casing.  The  input  and  output  shafts  are  only  coupled 
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(a) 


(b) 


Figure  4.16  Example  case  III:  (a)  Schematic  of  the  NASA  gear  test  facility.  (b)NASA 
spur  gear  pair  system  supported  on  rigidly  mounted  flexible  casing  [55]. 
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torsionally  to  the  rest  of  the  gear  test  facility  through  flexible  couplings,  and  the  geared 
system  is  driven  by  a  149kW  DC  motor.  Vibrational  level  at  various  positions  on  the 
casing  plate  are  measured  using  PCB  303 A  and  Endevco  2271  accelerometers  over  the 
operational  speed  range  of  2250-5750  rpm  which  corresponds  to  the  gear  mesh 
frequency  C0|i  range  of  1050-2683  Hz. 

4.7.2  Vibration  Models 

A  dynamic  finite  element  model  of  Figure  4.16b  is  developed  using  a  commercial 
software  [39]  as  outlined  in  Section  4.4.1.  The  magnitudes  of  k^y,  which  are  included 
in  the  FEM  model  through  a  6  dimensional  generalized  stiffness  matrix  element,  are 
computed  using  the  two  nonlinear  algebraic  equations  (4.27)  which  neglect  the  static 
elastic  deformation  of  the  casing  plate.  The  shafts  are  modeled  using  2  noded 
Timonshenko  beam  elements  with  axial  degrees  of  freedom  capability.  Four  noded 
quadrilateral  plate  elements  with  shear  deformation  and  rotary  inertia  effects  are  used  to 
construct  the  flexible  casing.  The  shafts  and  torsional  couplings  exterior  of  the  test 
gearbox  are  modeled  as  torsional  stiffness  elements  only,  and  the  gear,  pinion,  motor 
and  load  are  incorporated  in  the  FEM  model  using  generalized  mass  and  inertia 
elements.  Gear  mesh  coupling  between  the  gear  and  pinion  is  described  by  a  generic 
stiffness  matrix  of  dimension  6,  similar  to  Free  rotational  O^ft)  boundary 

condition  is  specified  at  the  motor  and  load  inertias,  while  ideally  clamped  uwa=Qwa=0, 
w=x,y,z,  boundary  points  at  the  comer  of  the  bottom  plate  are  assumed  to  simulate  the 
rigid  mount  feet. 

Natural  frequencies  and  modes  predicted  by  a  similar  FEM  model  have  already 
been  verified  by  comparison  with  experimental  modal  analysis  in  an  earlier  publication 
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by  us  [50].  Hence,  eigensolutions  will  not  be  discussed  here.  But,  it  may  be  noted  that 
the  conventional  geared  rotor  dynamic  models  with  simple  bearing  models  and  without 
casing  dynamics  are  not  valid  for  high  excitation  frequencies  (beyond  the  first  casing 
plate  elastic  mode).  For  the  forced  response  study,  only  the  fundamental  harmonic  ©h 
of  the  static  transmission  error  excitation  e(t)  =  3.5  sin(©{,t)  pm  computed  from  the  gear 
tooth  profile  and  geometry  [46].  About  110  dynamic  degrees  of  freedom  are  selected  to 
minimize  computational  effort  while  still  maintaining  sufficient  accuracy.  Modal 
damping  ratio  £  =  0.05  is  assumed  over  the  frequency  range  of  interest 


4.7.3  Casing  Response 

Several  locations  on  3  different  casing  plate  surfaces  as  shown  in  Figure  4.17  have 
been  chosen  for  the  experimental  validation  of  our  theory.  Figure  4.18  compares 
predicted  and  measured  mean  square  acceleration  spectra  at  ©j,  over  0  kHz  to  4  kHz 
even  though  measurements  have  been  conducted  only  between  1  kHz  to  3  kHz.  In 
general,  good  agreement  is  found  between  experiment  and  theory.  Some  discrepancy  is 
observed  above  2.4  kHz  which  is  due  to  the  limitation  of  the  FEM  model  in  the  high 
modal  density  regime.  It  may  be  noted  that  the  simple  theory  which  utilizes  only 
conventional  bearing  models  cannot  predict  flexural  vibration  of  the  casing  plate  as 
shown  in  Figure  4.18.  A  broad  band  vibratory  energy  comparison  is  achieved  by 
averaging  the  mean  square  value  of  the  acceleration  over  one-third  octave  bands,  in 
dB  per  unit  frequency  bandwidth  A©,  which  is  defined  as  follows 


L  ,dB  =  10  log 


10 


(isr)2> 

A  a) 


(aa*)  1 


<  a2> 


ref 


;  <A2>ref  =  1.0  g2  (4.35) 
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Figure  4. 17  Measurement  locations  on  the  casing  plate  for  example  case  HI. 


Acceleration  (g)  Acceleration  (g)  Acceleration  (g) 


149 
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Figure  4.18  Comparison  of  casing  flexural  vibrations  for  example  case  III  with  8zm  = 
0.08  mm.  (a)  Side  plate  supporting  the  bearings  and  facing  the 
dynamometer  (location  18a).  (b)  Top  plate  cover  (location  1 8b).  (c)Side 
plate  with  no  bearings  (location  18c). 


& 
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where  A*  implies  the  complex  conjugate  of  the  acceleration.  Our  predictions  are  within 
+5.0  dB  of  the  measured  values  as  shown  in  Table  4.8. 

The  sensitivity  of  results  to  [K]bm  is  evaluated  next.  Recall  that  [K]bm  is 
computed  based  on  the  knowledge  of  the  mean  axial  displacements  $zm  applied  to  the 
bearings.  This  parameter  has  been  assumed  to  be  a  constant  in  our  analysis  although  its 
exact  value  for  the  NASA  experimental  setup  is  not  known.  Using  2  reasonable  values 
of  8zm,  we  obtain  a  range  of  acceleration  spectra  in  Figure  4.19.  Almost  all  the 
experimental  data  are  now  within  this  range.  This  explains  experimental  scatter 
observed  in  Figure  4.18  as  experimental  S2in  may  vary  slightly  from  one  steady- state 
speed  to  another  during  testing. 

In  example  case  II,  we  found  that  only  a  few  coefficients  in  [K]bin  are  necessary 
depending  on  the  excitation.  Accordingly,  we  had  used  only  kbyy,  kby0x  and  kb0x0x. 
But  in  the  present  FEM  model  of  case  III,  we  must  include  the  entire  symmetric  [K]blii 
matrix  of  dimension  6  with  no  simplification  at  all.  Now  if  we  omit  the  off-diagonal 
term  kbx0y  which  is  not  directly  excited  by  e(t)  in  the  y-direction  of  Figure  4.1  while  still 

retaining  other  kbjj  coefficients,  no  significant  changes  in  our  predictions  are  seen. 
Similarly,  if  the  off-diagonal  terms  kbzy  and  kbz0x  which  are  related  to  the  shaft  axial 

degrees  of  freedom  are  also  neglected,  predicted  acceleration  spectra  vary  slightly. 
However,  if  kby0x  and  k^^  which  constraint  the  bending  motion  of  the  shaft  excited 

by  e(t)  are  assigned  zero  values  per  simple  theory,  zero  out-of-plane  acceleration  levels 
are  observed.  This  confirms  that  the  vibration  transfer  through  the  bearings,  from  the 
shaft  bending  motion  to  the  casing  flexural  motion,  is  highly  dependent  on  these  two 
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Table  4.8  Comparison  of  (dB)  for  example  case  HI. 


Center  Location  in  Figure  4. 1 7 

Frequency  (Hz)  * - * - 1 

18(a)  18(b)  18(c) 


1250 

Theory 

-3.0 

-13 

-4.5 

Experiment 

-2.0 

-17 

-8.4 

1600 

Theory 

-5.5 

-7.0 

-11 

Experiment 

-10 

-8.1 

-14 

2000 

Theory 

-11 

-5.4 

2.1 

Experiment 

-10 

-7.4 

0.0 

2500 

Theory 

-3.0 

0.5 

5.9 

Experiment 

-7.3 

-1.0 

4.0 

52 


Frequency  (Hz) 


F‘Surc  4-19  Sensitivity  of  predicted  casing  acceleration  to  preload  8^  for  example  case 
III.  (a)  Side  plate  supporting  the  bearings  and  facing  the  dynamometer 
(location  18a).  (b)  Top  plate  cover  (location  18b).  (c)  Side  plate  with  no 
bearings  (location  18c). 
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4.8  CONCLUDING  REMARKS 

A  generic  geared  rotor  system  model  has  been  developed  using  lumped  parameter 
and  dynamic  finite  element  techniques  which  incorporates  a  new  mathematical  model  for 
precision  rolling  element  bearings  proposed  earlier  in  Chapter  II.  This  system  oriented 
model  includes  internal  rotating  system,  rolling  element  bearings,  flexible  torsional 
couplings,  motor  and  load,  flexible  or  rigid  casing,  and  compliant  or  massive  mounts. 
The  discrete  shaft  model  excludes  the  effect  of  gyroscopic  moment,  but  includes 
torsional,  flexural  and  longitudinal  motions.  In  example  cases  I  and  II,  only  flexural 
and  torsional  motions  are  predicted.  The  effects  of  mount  stiffness,  casing  mass  and 
bearing  preload  on  the  overall  dynamic  behavior  have  been  investigated  through  3 
example  cases  of  single-stage  rotor  and  geared  rotor  systems  excited  by  rotating  mass 
unbalance  at  shaft  frequency  and  kinematic  transmission  error  excitation  at  mesh 
frequency  respectively.  The  results  indicate  that  our  proposed  model  is  more  accurate 
than  conventional  models  given  in  the  literature.  For  instance,  we  are  able  to 
predict  bearing  and  mount  moment  transmissibilities  and  improve  casing  flexural 
vibration  prediction  significantly  using  our  theory;  measurements  made  on  case  III 
validate  our  formulations  partially.  Conversely,  the  conventional  models  are  unable  to 
account  for  any  bearing  moment  transmissibility,  and  consequently  predict  zero  flexural 
response  on  the  casting. 

Natural  frequencies  of  the  example  case  III  geared  rotor  system  increase  with 
higher  mount  stiffness  and  bearing  preload,  and  lower  casing  mass.  Transition  mass 
and  stiffness  points  which  separate  regimes  of  a  natural  frequency  diagram  are 
predicted.  Through  extensive  parametric  studies,  we  find  that  bearing  transmissibilities 
may  be  reduced  by  using  a  higher  bearing  preload.  However,  this  may  result  in  higher 
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dynamic  transmission  error  for  the  spur  gear  pair,  mostly  due  to  the  coupled  flexural- 
torsional  motion.  Mount  transmissibility  is  affected  by  changes  in  casing  mass,  mount 
stiffness  and  mean  bearing  axial  displacement,  but  bearing  transmissibility  is  most 
sensitive  to  the  mean  bearing  axial  displacement.  Similar  conclusions  can  be  drawn 
regarding  the  other  two  example  cases. 

Our  theory,  though  restricted  to  the  linear  and  time-invariant  dynamic  system,  is 
comprehensive.  It  can  be  used  for  analysis  as  well  as  design  studies  of  other  rotating 
mechanical  systems  with  multiple  shafts  and  gear  pairs  or  multi-staged  rotors. 


CHAPTER  V 


STATISTICAL  ENERGY  ANALYSIS 

5.1  INTRODUCTION 

It  has  been  shown  in  Chapters  H-IV,  using  classical  lumped  parameter  and  dynamic 
finite  element  techniques,  that  the  proposed  bearing  model  is  more  accurate  than 
existing  models  for  predicting  vibration  transmission  through  bearings  in  a  geared  rotor 
system.  Although  the  proposed  model  has  been  shown  to  be  reliable  up  to  a  moderately 
high  frequency,  it  is  conceivable  that  this  model  is  inadequate  at  very 
high  frequencies  where  the  modal  density  is  high.  Classical  vibration  models  do  not 
predict  modes  accurately  in  this  freouency  regime,  and  even  if  it  is  possible  to  do  so  by 
employing  closely  spaced  nodal  points,  such  models  require  a  significantly  large 
computational  effort  Moreover,  the  vast  amount  of  predicted  response  spectra  at  many 
spatial  points  would  be  difficult  to  interpret.  Accordingly,  asymptotic  or  statistical 
methods  must  be  adopted;  typical  techniques  include  the  statistical  energy  analysis  [56- 
60],  asymptotic  modal  analysis  [61-63]  and  asymptotic  analysis  using  infinite  system 
impedances  [64]. 

This  study  concentrates  on  the  development  of  a  broad  band  vibratory  energy 
transfer  model  for  a  geared  rotor  system  with  the  proposed  bearing  model  using  the 
statistical  energy  analysis  (SEA)  method.  This  method  has  been  applied  successfully  to 
a  wide  variety  of  structural  dynamic  and  acoustic  systems  with  large  number  of  modes 
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[56-60,65-69].  However  several  unresolved  research  issues  still  exist  [57,61,62,70, 
71].  The  specific  objectives  of  this  study  are  to:  (i)  conduct  modal  analysis  of  a  geared 
rotor  system  with  flexible  casing  and  mounts,  (ii)  investigate  the  feasibility  of  applying 
SEA  to  this  problem,  (iii)  analyze  the  following  4  example  cases  using  SEA:  I.  a  plate- 
cantilevered  rectangular  beam,  n.  case  I  with  circular  shaft-bearing  system  replacing  the 
cantilevered  beam,  HI.  a  circular  shaft-bearing-plate-mount  system,  and  IV.  a  simple 
geared  rotor  system,  and  (iv)  perform  parametric  studies  to  examine  the  characteristics 
of  vibratory  energy  transfer  through  bearings,  and  mean  square  vibroacoustic  response 
of  the  casing.  The  first  and  second  example  cases  are  revised  and  extended  versions  of 
a  study  performed  by  Lyon  and  Eichler  [56,58].  The  later  two  example  cases  are  taken 
from  Chapters  III  and  IV  where  these  systems  have  been  studied  at  low  frequencies. 
Experimental  validations  are  also  included. 

5.2  MODAL  ANALYSIS  OF  GEAR  CASING  AND  MOUNTS 

Analytical  and  experimental  modal  analyses  of  a  real  gear  casing  and  mounts  have 
been  conducted  to  investigate  the  feasibility  of  using  SEA.  Natural  frequencies  coj  and 
modes  4>j  are  calculated  using  a  commercial  finite  element  method  (FEM)  program  [39], 
and  predictions  for  a  rigidly  mounted,  stiffened  gearbox  are  verified  by  the  experimental 
modal  analysis  (EMA). 

5.2.1  Finite  Element  Model 

The  rectangular  gearbox  as  shown  in  Figure  5.1a  is  approximately  0.254m  x 
0.279m  x  0.330m  (lO.Oin  x  1 1 ,0in  x  13.0in),  and  all  of  its  plates  are  0.006m  (1/4  in) 
thick  made  of  1020  steel  excej  t  the  regions  near  the  bearings  which  are  0.025m  (l.Oin) 


INPUT  SHAFT 


Flexible 

mount 


0.0064m  thick 
aluminium  fuselage 
sheet  with  dimensions 
0.7620m  x  0.6640m 
supported  in  between 
gearbox  and  mounts 


(b) 


Figure  5. 1  Gearbox  with  rigid  and  flexible  mounts. 
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thick.  There  are  four  circular  holes  for  the  bearings,  two  at  each  side  plate  supporting 
the  shafts.  Figure  5.1b  illustrates  the  0.254m  (lO.Oin)  tall  flexible  mount  frame  which 
is  constructed  from  eight  0.006m  (1/4  in)  thick,  1020  steel  angle  beams  with  three 
different  lengths.  Fuselage  sheet,  as  shown  in  Figure  5.1b,  is  attached  horizontally  to 
the  flexible  mount  structure.  The  casing  is  supported  at  each  comer  of  the  base  plate  for 
all  mounting  conditions,  and  the  mounts  are  attached  to  a  rigid  foundation. 

Two  FEM  models  of  the  rectangular  gearbox  without  its  spur  gears  set,  shafts  and 
bearings  are  developed  for  the  rigidly  and  flexibly  mounted  casings.  These  FEM 
models  consist  of  four-noded  quadrilateral  plate  elements  with  bending  and  membrane 
capabilities  for  the  housing  and  attached  fuselage,  and  two-noded  shaft  element  with 
shear  deformation  and  rotary  inertia  capabilities  for  the  flexible  mount  skeleton  and 
housing  plate  stiffeners.  The  boundary  conditions  are:  (i)  zero  displacements  and 
rotations  at  each  comer  of  the  base  plate  for  the  rigid  mount,  and  (ii)  similar  conditions 
at  each  foot  of  the  flexible  mount.  The  interfaces  between  adjacent  housing  plates  are 
assumed  to  be  continuous.  About  100  dynamic  degrees  of  freedom  are  specified  to 
reduce  computational  effort  while  still  maintaining  sufficient  accuracy.  Natural 
frequencies  are  computed  up  to  at  least  4  kHz  to  cover  the  gear  mesh  frequency  regime. 

5.2.2  Experiments  and  Model  Validation 

Modal  experiments  have  been  performed  on  a  NASA  high  precision  gearbox  with 
the  spur  gear  set,  shafts,  and  ball  bearings  installed.  An  approximate  configuration  of 
the  NASA  gearbox  is  shown  in  Figure  5.2.  The  nominal  dimensions  of  the  gear 
housing  have  been  given  in  the  previous  section.  The  variable  center  distance  gear-shaft 
pair  is  supported  by  four  ball  bearings.  Four  side  plates  and  a  base  plate  are  welded 
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together  while  the  top  plate  is  bolted  to  the  side  plates.  The  housing  plates  are  stiffened 
internally,  and  the  gear  housing  system  is  mounted  rigidly  to  a  massive  foundation. 
Dynamic  transfer  functions  are  obtained  only  on  the  exterior  of  the  gear  casing  structure 
using  the  GENRAD  2515  system  [72].  For  these  experiments,  154  measurement  points 
have  been  selected  in  the  direction  transverse  to  the  plane  of  the  casing  plates  with  the 
reference  point  being  approximately  near  the  center  of  the  top  plate  to  avoid  nodal  points 
of  interest.  Natural  frequencies  and  modes  are  then  estimated  using  the  MODAL  PLUS 
program  [73].  Here,  the  exponential  method  has  been  used  to  extract  modal  parameters 
and  generate  analytical  functions  for  the  transfer  functions,  while  the  circle  fit  method 
has  been  used  to  construct  the  modal  vectors. 

Figure  5.3  compares  predicted  and  measured  coj,  and  EEM  is  found  to  be  in  good 
agreement  with  EMA.  For  each  mode  <|>j,  two  simplified  illustrations  are  shown  in 
Figure  5.3:  (i)  mode  shape  of  the  three  visible  plates  .and  (ii)  mode  shape  of  the  three 
nonvisible  plates  in  an  approximate  isometric  view.  The  higher  modes,  not  shown  here, 
are  also  given  by  similar  combinations  of  plate  flexural  motions.  Comparison  between 
theory  and  experiment  for  the  higher  modes  is  made  on  the  basis  of  number  of  modes 
within  each  one-third  octave  bands  over  500-4000  Hz  in  Figure  5.4a,  because  of  the 
high  number  of  participating  modes  observed.  The  results  again  indicate  that  FEM  is  in 
good  agreement  with  EMA. 

5.2.3  Parametric  Studies 

The  rigidly  mounted  gear  housing  is  observed  to  possess  only  the  elastic  modes  of 
the  casing  plates.  On  the  other  hand,  FEM  model  of  the  flexibly  mounted  gear  casing 
indicates  that  the  first  six  modes  are  translational  and  rotational  rigid  body  modes  of  the 
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(0,  =  501  Hz  co 2  =  598  Hz 
(497  Hz)  (584  Hz) 


co 3  =  627  Hz  co4  =  752  Hz 

(639  Hz)  (781  Hz) 


co  5  =  889  Hz 
(877  Hz) 


Figure  5.3  Mode  shapes  of  the  rigidly  mounted,  stiffened  NASA  gearbox  as  predicted 
by  FEM  (EMA  results  are  given  in  parenthesis).  Bold  sign  implies  larger 
amplitude. 
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One-Third  Octave  Band  Center  Frequency  (kHz) 


(a) 


<b) 


Figure  5.4  Number  of  modes  within  each  one-third  octave  bands  for  the  NASA 

gearbox,  (a)  Predicted  and  measured  number  of  modes  for  rigidly  mounted 
casing,  (b)  Number  of  modes  for  rigid  and  flexible  mounting  conditions. 
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the  casing  as  shown  in  Figure  5.5;  for  example,  £0j  =  54  Hz  corresponds  to  the  casing 
vibrating  in  the  Y-direction  as  shown.  In  addition,  (Oj  are  considerably  lower,  by 
approximately  one  order  of  magnitude  as  compared  to  those  of  the  rigidly  mounted  gear 
casing.  These  rigid  body  vibration  modes  result  from  the  complex  elastic  deformations 
of  the  flexible  mount  skeleton  and  fuselage  sheet.  The  casing  plate  natural  frequencies 
are  also  lowered,  especially  the  first  few,  when  the  box  is  mounted  flexibly.  Figure 
5.4b  compares  the  number  of  modes  of  the  flexibly  mounted  gear  casing  to  that  of  the 
rigidly  mounted  one.  High  modal  counts  above  1  kHz  are  seen. 

The  introduction  of  gear  casing  plate  stiffeners  as  shown  in  Figure  5.2  does  not 
change  the  nature  of  the  mode  predictions.  Also,  the  natural  frequencies  for  this  case 
only  vary  slightly;  lower  modes  are  affected  more  by  the  stiffeners  than  the  higher  ones. 
Note  that  the  numbers  of  modes  in  each  1/3  octave  band  over  400-4000  Hz  range  remain 
nearly  the  same  as  evident  from  Table  5.1. 


Table  5. 1  Number  of  modes  for  the  stiffened  and  unstiffened  gearboxes  mounted 
rigidly. 


1/3  Octave  Band  Center 

Frequency  (kHz) 
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Figure  5.5  Rigid  body  modes  of  the  flexibly  mounted  gearbox  as  predicted  by  FEM. 
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5.3  JUSTIFICATION  FOR  USING  SEA 
5.3.1  Modal  Densities 

Analytical  and  experimental  modal  analyses  of  the  NASA  gear  casing  clearly  show 
that  modal  density  tends  to  be  high  at  higher  frequencies.  Although  the  analysis  is 
restricted  to  a  gear  casing  system,  it  is  reasonable  to  assume  that  similar  results  are  valid 
for  other  rotating  mechanical  system  housings.  Accordingly,  SEA  can  be  used 
justifiably. 

Next,  consider  the  modal  densities  of  a  shaft  or  beam  ns  and  rectangular  plate  nc, 
given  for  bending  motion  with  simply-supported  boundary  conditions  [56] 


/jv iti 

4/  4  EI  s  to2 


Ac  / 3pc  (1  ~  M-2) 
nc“  hcV  E 


(5-1) 


where  p  is  the  material  density,  Ac  is  the  plate  surface  area,  Ls  is  the  shaft  length,  E  is 
the  modulus  of  elasticity,  Is  is  the  area  moment  of  inertia  of  the  shaft,  co  is  the 
bandwidth  center  frequency,  h  is  the  plate  or  shaft  thickness,  p  is  the  Poisson's  ratio, 
and  the  subscripts  s  and  c  denote  shaft  and  plate  respectively.  For  typical  numerical 
values  given  later  in  Section  5.4  we  find  that:  (i)  nc  =  0.227,  and  (ii)  ns  =  0.012  at  100 
Hz,  ns  =  0.004  at  1000  Hz  and  ns  =  0.002  at  4000  Hz.  Note  that  ^  »  ns  as  expected. 
Although  the  shaft  has  a  fairly  low  number  of  modes  in  the  frequency  range  of  interest, 
SEA  is  still  valid  since  the  plate  modal  density  nc  is  very  high. 
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5.3.2  Literature  Review 

A  majority  of  publications  are  on  the  application  of  the  SEA  procedure  to 
dynamic  systems  with  high  n  such  as  structural-acoustical  interactions  in  a  fuselage 
[56,65,74],  sound  transmission  through  panels  [60],  and  vibratory  energy  transmission 
in  mechanical  equipment  [56].  Of  interest  here  are  analytical  or  experimental  estimation 
of  SEA  parameters  for  simple  structural  systems  described  by  a  flat  plate,  a  cylinder 
and/or  a  shaft  [56,58,59,66].  In  these  studies,  structural  connections  are  often  assumed 
rigid  such  as  in  the  ideally  welded  case. 

Lyon  and  Eichler  [56,58]  in  1964  and  Lyon  and  Scharton  [59]  in  1965  developed 
analytical  expressions  for  the  coupling  loss  factor  rj  in  several  connected  structures, 
such  as  a  plate  bonded  to  a  cantilevered  beam  [56,58].  Here  tj  was  derived  assuming  a 
semi-infinite  beam  attached  to  an  infinite  plate  and  by  further  assuming  that  only  a 
dynamic  moment  coupling  at  the  joint  can  describe  the  motion/force  transmission 
phenomenon.  This  problem  is  re-examined  in  this  study  and  is  then  extended  to  a 
circular  shaft-bearing-plate-mount  system.  In  addition,  Lyon  and  Eichler  [56,58]  also 
developed  SEA  models  of  two  structures  inter-connected  through  a  single  (scalar) 
stiffness  element.  Typical  examples  include  two  longitudinal  rods  connected  via  a  linear 
stiffness  element,  two  discrete  masses  coupled  by  a  linear  spring,  and  a  plate  attached  to 
a  single  degree  of  freedom  resonator  [56].  A  similar  analysis  on  the  longitudinal 
vibration  of  linearly  coupled  rods  was  performed  by  Keane  and  Price  [69]  in  1987. 
Loss  factors  of  typical  line  or  point  connected  structures  such  as  a  plate  welded  to  a 
cylinder,  cross  beams  and  two  perpendicular  plates  bolt  .u  or  welded  together,  have  been 
calculated  assuming  ideal  rigid  joints  [56,58,59,75,76].  But  a  compliant  bearing  system 
problem  is  yet  to  be  analyzed. 
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5.4  EXAMPLE  CASE  I:  COUPLING  LOSS  FACTOR  OF  PLATE- 
CANTILEVERED  BEAM  SYSTEM 

First,  we  attempt  to  rework  the  plate-cantilevered  beam  problem  of  Lyon  and 
Eichler  [56,58]  as  shown  in  Figure  5.6.  Only  flexural  motions  of  the  plate  and  beam  are 
considered  in  this  case.  Accordingly,  Lyon  and  Eichler  [56,58]  developed  an 
expression  for  the  coupling  loss  factor  r\sc  which  describes  the  vibratory  energy  transfer 
between  the  beam(s)  and  the  plate(c)  due  to  a  moment  coupling  at  the  joint 


(2PsAskscs) 
torn  „ 


Re 


(5.2) 


where  k=a/vas 's  the  radius  of  gyration,  c=V  E/ps  is  the  wave  speed,  m  is  the  mass 
and  Re(  )  implies  the  real  part  of  a  complex  variable.  The  driving  point  moment 
impedances  for  the  plate  Zc  and  beam  Zg  are  [58,77] 


Z  c((o)  = 


16pchciC  c 


-  ,  Ps  AsKscs^s  ,  .  .  . 

.  Zs((0)  = - - - (1-1) 


(5.3) 


where  k  is  the  wavenumber  and  h  is  the  plate  or  beam  thickness.  Here,  note  that  Zs  of 
equation  (5.3)  differs  from  Lyon  and  Eichler's  expression  [56,58]  in  the  sign  of  the 
imaginary  part  which  is  probably  a  typographical  error.  They  assumed  that  IZCI»IZSI, 
ks=kc,  Cs=cc  Ps=Pc.  and  simplified  equation  (5.2)  to  yield  a  frequency  invariant 
expression  for  rise  as 


t1sc  =  Ws/(4Ls) 


(5.4) 
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Figure  5.6  Schematic  of  example  case  I:  a  plate-cantilevered  beam,  (a)  Experimental 

finite  System,  (b)  Theoretical  infinite  system.  i 

t 

v 
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where  Ws  and  Ls  are  the  beam  width  and  length  respectively  as  shown  in  Figure  5.6. 
However,  our  calculations  indicate  that  IZCI>IZSI  but  not  IZCI»IZSI  for  the  parameters 
used  by  Lyon  and  Eichler  [58],  as  shown  in  Figure  5.7a.  Using  equation  (5.2),  T|sc  is 
recomputed  and  compared  with  the  approximate  model  given  by  equation  (5.3)  in  Figure 
5.7b.  Experimental  results  given  by  Lyon  and  Eichler  are  also  plotted  in  Figure  5.7b.  It 
can  be  seen  that  our  calculation  is  better  than  earlier  prediction.  At  low  frequency,  a 
discrepancy  is  observed  between  theory  and  experiment  which  is  primarily  due  to  the 
low  modal  density  in  this  regime.  The  presence  of  a  low  natural  frequency  may  be  due 
to  the  compliant  epoxy  bond  between  the  beam  and  plate.  However,  above  the 
threshold  frequency  where  many  modes  participate,  shown  as  a  vertical  line  in  Figure 
5.7b,  the  slope  of  the  least  square  straight  line  fit  on  the  experimental  data  is  nearly  the 
same  as  the  predicted  r]^. 

5.5  EXAMPLE  CASE  II:  COUPLING  LOSS  FACTOR  OF  CIRCULAR 

SHAFT-BEARING-PLATE  SYSTEM 

Next,  we  modify  Figure  5.6  by  inserting  a  ball  bearing  between  the  circular  shaft 
(which  replaces  the  beam  in  Figure  5.6)  and  the  rectangular  plate.  Again,  a  semi-infinite 
shaft  and  an  infinite  plate  are  assumed.  For  SEA,  we  reduce  the  system  to  a  plate 
subsystem  and  a  shaft-bearing  subsystem.  The  coupling  loss  factor  T)sc  is  still  given  by 
equation  (5.2),  but  Zj  must  be  modified  to  account  for  the  compliant  bearing. 

Consider  a  shaft  with  boundary  conditions  shown  in  Figure  5.8.  The  bearing  end 
is  subjected  to  zero  transverse  velocity  vysa(t,z=0')=0  and  a  sinusoidally  varying 
moment  Mxsa(t,z=0')=MxsaeioM.  Here,  z=0'  refers  to  the  junction  point  between  bearing 
and  plate,  and  z=0+  is  the  junction  point  between  the  bearing  and  shaft.  In  the 
frequency  domain,  using  the  definition  of  bearing  force  vector 


Coupling  Loss  Factor, 
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Figure  5.7  Comparison  between  Lyon  and  Eichler's  approximation  [56,58]  and  our 
proposed  formulation  for  example  case  I.  (a)  Comparison  of  \Z^(\ZC+Z^. 
(b)  Comparison  of  predicted  rjsc  with  experimental  results  given  by  Lyon 
and  Eichler  [56,58]. 


171 


Mxsa(W>Z=0  )  =  Mxsa 
vysa(w,  z=0")=° 


z=0 


Figure  5.8  Boundary  conditions  for  example  case  II:  semi-infinite  shaft-bearing 
system. 
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{f(CO)}b=[K]b  m{v(o))}b/io)  in  terms  of  the  velocity  vector  {v(co)}b  and  the  proposed 
bearing  stiffness  matrix  [K]bm,  vysa(co,z=0‘)=0  is  shown  to  be  equivalent  to 
Eysa(w>0  )=(^byex/^b9xex)^xsa+vysa(°i>0+)  {  (^by0xAb6x0x^^bex0x“^byy  }  / 0W)- 
The  governing  equations  for  Mxsa(co,z)  and  Fysa(co,z)  of  the  shaft  in  terms  of  vysa(co,z) 
are 


M 


xsa 


El  s  ^  vysa 
dz  ^ 


El  s  d  vysa 
l05  dz3 


(5.5) 


In  general,  the  bounded  solution  for  vysa(co,z)  is 


vysa(co,z)  =  {  bj  e  +  b2  e  k*z  } 


(5.6) 


Using  equations  (5.5)  and  (5.6),  we  enforce  the  boundary  conditions  for  the  shaft  at 
z=0+  to  obtain  the  following  closed  form  solution  for  the  coefficient  vector 
(b)  =  {b1,b2}T. 
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Hence,  the  bearing-shaft  impedance  Zg  is  then  given  by 


Z  c(co)  = 


MXSa(W,°) 


bexex 


dv 


ysa 


dz 


(£0,0") 


ICO  + 


k.  Q  v  (£o,0+)+k 
by0  ysa 


dv 


(5.8) 


1)0,0,  dz 


ySa  (to,  0+) 


M 


xsa 


Although  Mxsa  appears  in  equation  (5.8),  Z ^  is  independent  of  the  magnitude  of  Mxsa 
dV  . 

since  vysa(co,0+)  and  — ^ — (co,  0  )  are  linearly  proportional  to  in  equation  (5.5). 

Using  equations  (5.2),  (5.7)  and  (5.8),  T]sc  is  computed  for  several  bearings 
whose  stiffness  coefficients  are  tabulated  in  Table  5.2.  Figure  5.9  compares  these 
where  the  material  and  geometrical  properties  of  the  shaft  and  plate  are  equivalent  to 
those  used  by  Lyon  and  Eichler  [56,58]  in  example  case  I.  For  a  very  soft  bearing,  Zs 
is  dominated  by  the  bearing  parameters,  and  hence  \7^\»\Z^  is  valid.  This  implies  that 
T]x  reduces  to  the  frequency  invariant  expression  given  by  equation  (5.4)  as  shown  by 
set  A  in  Figure  5.9.  On  the  other  hand,  risc  for  a  rigid  bearing  (say  set  B  or  C)  is 
typically  smaller  than  for  a  soft  bearing  given  by  A.  The  extreme  values  of  rise,  as 
kbij  or  bearing  preload  becomes  very  large,  depend  on  the  relative  magnitudes  of  the 
bearing  stiffness  coefficients  as  shown  in  Figure  5.9  for  sets  D  and  E.  Such  deviation 
for  a  stiff  bearing  from  set  B  is  primarily  due  to  the  presence  of  off-diagonal  stiffness 
coefficient  such  as  kjjye,.  However  it  is  found  that  for  most  stiff  bearings  will 
approach  set  E  in  Figure  5.9  which  is  obviously  identical  to  the  prediction  for  the 
example  case  I  using  equations  (5.2)  and  (5.3). 
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Table  5.2  Typical  bearing  stiffness  coefficients  of  Figure  5.9  for  example  case  II. 


Set 

kbyy  (N/m) 

^byG* (N) 

kb8x9x  (Nm) 

A 

Very  Compliant  Bearings  (kw:  <  1E3) 

B 

1E5 

5E4 

2E4 

C 

5E5 

8E4 

5E4 

D 

1E8 

1E7 

1E6 

E 

1E8 

3E7 

1E7 

5.6  EXAMPLE  CASE  III:  A  CIRCULAR  SHAFT-BEARING-PLATE 

SYSTEM 

5.6.1  Theory 

Consider  a  circular  shaft-bearing-plate-mount  system  similar  to  example  case  II  but 
with  a  shaft  of  finite  length  as  shown  in  Figure  5.10a.  Note  that  this  system  has  been 
analyzed  earlier  using  the  deterministic  vibration  modeling  technique  for  low  frequency 
response.  Recall  the  unconstrained  end  of  the  non-rotating  circular  shaft  is  subjected  to 
a  harmonically  varying  force  Fys(t)  =  Fysaeio)t  +  Fysm,  where  Fysm  is  the  mean  force. 
From  Chapter  III,  we  know  that  the  longitudinal  and  torsional  motions  of  the  shaft,  and 
the  in-plane  vibration  of  the  plate  can  be  neglected.  Hence,  the  shaft  bending  vibration 
and  plate  flexural  motion  are  of  interest  here.  Also,  the  previous  deterministic  vibration 
models  indicate  that  the  coupling  between  these  two  motions  is  mainly  due  to  the 
dynamic  moment  at  the  bearing  provided  the  longitudinal  shaft  vibration  is  not  excited. 
Accordingly,  two  subsystems  which  can  be  easily  identified  using  SEA  are  the 
transverse  modes  of  the  shaft-bearing  system  and  the  flexural  modes  of  the  plate-mount 
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Shaft-bearing  system 
(bending  modes) 
Es(to),  ns(to) 
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- ►  Energy  Dissipation 
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Energy  Transfer 
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Plate-mount  system 
(flexural  modes) 

Ec  (to),  nc  (to) 

Ac 

- ►  Energy  Dissipation 

(b) 


Figure  5.10 


Example  case  in:  (a)  schematic  of  the  circular  shaft-bearinc-plate-mouni 
system  and  (b)  A  SEA  model  of  the  shaft-bearing-plate-mlum  sysS 
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system  as  shown  in  Figure  5.10b.  We  now  proceed  to  derive  T|sc  governing  the 
vibratory  energy  flow  between  these  two  subsystems  using  the  driving  point  junction 
impedance  method  suggested  by  Lyon  [56]. 

Since  the  rectangular  plate  dimension  is  much  larger  than  the  bearing  dimension, 
the  plate  is  again  assumed  infinite.  Consider  the  vibrational  power  flow  nsc  from  the 
shaft-bearing  subsystem  (subscript  s)  to  the  plate  (subscript  c),  due  to  the  uniform 
external  Gaussian  random  force  over  a  frequency  bandwidth  Ao)  with  center  frequency 
to. 


11^(0))  =gjtisc(g))  ns(to) 


'  Es(to) 
,  ns(“) 


Ec(to)  ] 
nc(co)  J 


(5.9) 


where  nj  and  Ej  (j  =  s,  c)  refer  to  subsystem  modal  density  and  total  vibratory  energy 
respectively.  Since  the  plate  is  assumed  to  be  reasonably  well  damped  and  geometrically 
large,  equation  (5.9)  is  approximated  assuming  nc»ns  or  Ec/nc«Es/ns  to  yield 


nsc(ta) 
T1*(t0)  *  to  Es(to) 


(5.10) 


For  the  shaft,  Es  =  ms<Vs>  where  ms  is  the  shaft  mass  and  <Vg>  is  mean  square  shaft 

transverse  velocity  averaged  over  Aco  and  shaft  length  Lg.  Using  the  expression  for  11^ 
given  by  Lyon  and  Eichler  [58]  as  discussed  in  the  previous  section,  in  terms  of  Zg 
and  Zc  is  identical  to  equation  (5.2)  which  is  rewritten  here  for  the  circular  shaft  case  as 


Tlsc(<o) 


(5.11) 
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For  an  infinite  steel  plate  of  thickness  hc,  the  point  moment  impedance  Zc  is  given  by 
equation  (5.3);  it  is  rewritten  as 


'  ;  Ms<<1  <512> 


The  shaft-bearing  impedance  Zj  in  equation  (5.1 1)  is  derived  next  by  solving  the 
boundary-value  problem  for  the  system  shown  in  Figure  5.11a.  The  free  end  of  the 
shaft  has  a  vanishing  dynamic  shear  force  Fysa(co,z=Ls)=0  and  bending  moment 
MXSa(®,z=Ls)=0.  The  bearing  end  is  similar  to  Figure  5.8.  Following  the  same 

argument  used  previously  for  example  case  II  will  lead  to  the  two  boundary  conditions 
described  by  Mxsa(co,z=0+)=Mxsa  and  Fxsa(©,z=0+)  =  (kbx0y/kb0y0y)Mxsa  + 
vysa(©,O+){(kbx0y/kb0y0y)2kb0y0y-kbxx}/(i©).  Governing  equations  are  still  given 

by  equation  (5.5),  but  the  solution  for  vysa(co,z)  is  assumed  to  be  of  the  following  form 
[77] 


vysa(t,z)  =  {  &i  e  jk*z  +  t>2  e  jksZ  +  b3  e  k*z  +  b4  e  k*z  }  (5. 1 3) 

Using  equations  (5.5)  and  (5.13),  the  boundary  conditions  yield  the  following  algebraic 
problem 


(B]{b)  =  Mxsa{kby0x^cb0x0x,  io>/(EIkj),  0,  0) T 


(5.14) 


179 


^xsa(®’  z— 0  )  —  Mjjjj 

vysa(<o,z=(T)=0 


Mxsa(co,z=Ls)  =° 
Fysa(a),z=Ls)  =0 


0L 


0 - 

Bearing 


Flexible  Shaft 


-►  z 


z=0 


(a) 


vysa(«J,z=-L~)=0 

(£0,z=-L;)  =0 
dz  _ 

|  Bearing 


z=-Ls 


Mxsa(w,z=0  )=0 
Fysa(<o,z  =  0)  =Fxsa 


Flexible  Shaft 


(b) 


Figure  5.1 1  Boundary  conditions  for  example  case  HI:  finite  shaft-bearing  system. 

(a)  Moment  applied  at  the  bearing  end.  (b)  Force  applied  at  the  free  end. 
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The  nonzero  elements  of  coefficient  matrix  [B]  of  dimension  4  are 


Bll  =  £ 


B13=£ 


ikK  a 

3  byex 

Elk  :  -  ik .  +  1 

s  byy  k 


V 


bM*  j 


'  B12  =  ¥ 


ikh  n 

-  Elk  l  -  ik  +  -r— — 
s  byy  k 


pji  3  ,  A  by6x 

Elk  „  -  k,  + 


s  byy  k 


;  BU  =  ^ 


V 


-  Elk  „  -  kt  + 


bexexy 


bye, 


s  byy  k 


b®x®x  j 


®21~ ®22~ — ®23— — ®24-^  >  B-jj— — iB4j— ie  ^s^"s  ;  B^- i^42— — » 

B33=-B43=-e_ksLs  ;  B34=B44=eksLs  (5.15) 


Both  [B]  and  (b)  =  {bj,  b2,  b3,  b4}T  can  be  easily  obtained  numerically.  The  bearing- 

shaft  impedance  is  still  given  by  equation  (5.8). 

The  same  procedure  may  be  applied  to  obtain  the  driving  point  force  impedance  for 

a  harmonically  varying  transverse  force  Fysa(t,z=0).  Note  that  the  origin  is  redefined  at 

the  forcing  point  as  shown  in  Figure  5.1  lb  for  convenience.  The  boundary  conditions 

dv  vsa 

are  FySjj(o),z=0)— Fys3>  h4xg2(co,z=0)=0,  Vyjj(to,z=-Lj)-0,  and  ^  (to,  — Lj)=0. 
These  conditions  at  z=-L‘  can  be  rewritten  for  z=-L*  like  the  previous  case  to  obtain 

Mxsa((o,z:=— Lg)=— (kbQx0x  jT~(<a>  —  Es )  +kby0xvySa(co,— Ls ))  j (ico)  ; 

dv  w 

Bysa(®>z~— Ls)=— (kby0x  ^  (to,  —  Ls)  +  kbyyVySa(o),— Ls))^(ico)  (5.16) 

These  prescribed  boundary  conditions  again  yield  a  set  of  algebraic  problem  similar  to 
Equation  (5.14).  The  nonzero  elements  of  the  coefficient  matrix  [B]  of  dimension  4  are 
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Bn  -  -B21  -  B12  -  B22  -  -B13  -  iB23  -  -Bj4  -  -iB^  -  -1.0  ; 

B31  =  (“  jkskb8x0x  +  kby0)l  +  ^Dc  s  ^  e  > 

B32  =  ^kskb0x0x  +  kby0x  +  EDc  s  ^  e  ’ 

B 33  =  (  —  kskb0x0x  +  kby0x~  EDc  s  ^  e  ; 


B34  =  <kskb0x0x  +  kby0x-EIks)e"k‘L‘ 
B41  =  ^s^byfl!l  —  kbyy  —  s  ^e^C,L*  ; 
'42-'  “S,vby0x  “byy 


B,,  =  (-  Ut,kk..a  -kt„„  +  iEIk^)e  ik,L' 


-kb„  +  Elk3s>ek,L*  : 


B44  =  <-kSkbyerk»w'EIk“,e’k'L' 


(5.17) 


The  right  hand  side  vector  {b}  of  the  algebraic  problem  is  (0,Fysaco/(EIk^),0,0}T. 
Force  impedance  at  the  driving  point  is  then  given  by  Zs(co,z=0)=Fysa/vysa(co,0). 
Accordingly,  the  input  power  is  ns  =  (1/2)  F^  Re{(l/Zg)*}  where  Re{  }  is  the  real 
part  of  the  complex  variable  and  ( )*  implies  the  complex  conjugation. 

We  can  now  compute  the  vibratory  energy  transfer  11^  through  the  bearing  and 
steady-state  subsystem  energy  levels  Es  and  Ec  by  applying  the  energy  balances  to  both 
subsystems  show’n  in  Figure  5.10b;  here  Tigs  =  nji^. 


't\s+t\x((0) 

-lies  (C0) 

fns/“. 

(5.18) 

-Ilscfto) 

llc  +  W®). 

lEc(co)J  1 

1  o  . 
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Es«d) 

Ec(C0) 


_ ns  ( c  ^  cs) _ 

“(^c  +  ^sc  +  ^cs) 

_ ns  sc _ 

“(ns^c+Tic^sc+nsncs) 


(5.19) 


Since  E  =  m  <V2>,  the  following  velocity  levels  may  be  obtained  at  any  center 
frequency  (0  from  either  equation  (5.18)  or  (5. 19' 

<  Vc  _  m  s  (  21  (S  201 

<Vs2>‘mcUc+^(«»J 

5.6.2  Validation  and  Parametric  Studies 

In  order  to  validate  our  SEA  formulation,  we  compare  the  mean  square  mobility 

level  of  the  plate  with  experimental  data  reported  earlier  in  Chapter  III.  Note  that 

although  all  nonzero  bearing  stiffness  coefficients  kby  are  computed  and  given  in  Table 

3.6,  only  kbyy=3.69E8,  kbygx=3.52E5  and  kbgx0x=4.19E4  are  used  as  they  appear  to 

be  the  most  significant  ones  according  to  the  proposed  theory.  Using  equation  (5.19), 
<Vc>  is  computed  and  compared  with  experimental  results  in  Figure  5.12.  Theoretical 

predictions  for  three  values  of  dissipation  loss  factor  rjs=r|c5itrjs(co)  are  given  since  the 
choice  of  structural  damping  is  critical  to  the  SEA  analysis.  It  can  be  seen  from  Figure 
5.13  that  the  experimental  data  are  approximately  bounded  by  rjs=0.0003  and  t|s=0.03. 
Here  t|s  is  also  assumed  to  be  frequency  invariant  which  may  not  be  valid  in  the 
experiment.  Accordingly,  comparison  between  theory  and  experiment  is  deemed  to  be 
excellent. 

Further  comparison  between  theory  and  experiment  can  be  made  for  the  case  of  a 
semi-infinite  shaft  considered  in  the  example  case  II.  Using  equations  (5.2),  (5.7)  and 


1/3  Octave  Band  Center  Frequency  (Hz) 


Figure  5.12  Comparison  between  theory  and  experiment  [17]  foi  example  case  in  with 

2 

a  finite  shaft.  Here  plate  mobility  level  in  dB  is  101og|o<Vc>  re 
<V?>=1.0(m/s)2  for  FV„9=1.0N. 
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(5.8)  for  rise,  the  mobility  levels  are  computed  and  are  found  to  be  given  by  straight 
lines  as  shown  in  Figure  5.13.  These  lines  represent  the  asymptotic  behavior  of  the 
system  when  the  shaft  is  very  long  i.e.  Ls  — >  °o.  Also,  note  that  Figure  5.13  is 
consistent  with  the  trends  reported  by  Lyon  and  Eichler  [58]  for  the  plate-cantilevered 
beam  problem.  Again,  most  experimental  data  are  bounded  within  the  range  given  by 
ris  =  0.0003  to  0.03. 

Consider  the  finite  shaft  length  1^=1. 32  m  of  high  modal  density  ns.  In  Figure 
5.14  ri^  is  compared  for  a  finite  shaft  and  semi-infinite  shaft  It  can  be  seen  that  the 
result  for  the  semi-infinite  shaft  follows  the  average  values  of  the  finite  shaft. 

Next,  the  effect  of  bearing  preload  or  stiffness  coefficients  on  the  mean-square 
velocity  response  of  the  plate-mount  system  is  evaluated.  Equation  (5.20)  is  used  to 
predict  the  curves  shown  in  Figure  5.15  while  keeping  other  system  parameters  the 
same.  The  bearing  coefficients  used  in  this  analysis  arc  from  Tables  3.3  and  3.4.  Three 
bearing  preloads  used  here  are  115N,  190N  and  285N  which  are  referred  to  as  low, 
medium  and  high  preloads.  We  observe  minor  changes  with  preload  except  in  the 


,  vicinity  of  4  kHz. 


100  1000  10000 

1/3  Octave  Band  Center  Frequency  (Hz) 


Figure  5.15  Effect  of  bearing  preload  on  the  normalized  casing  plate  response  (example 
case  III). 
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5.7  EXAMPLE  CASE  IV:  A  GEARED  ROTOR  SYSTEM 
5.7.1  Assumptions 

As  the  final  case,  we  investigate  a  generic  geared  rotor  system  with  4  bearings  as 
shown  in  Figure  5.16a.  This  system  is  not  only  a  natural  extension  to  example  in,  it 
can  also  be  treated  as  an  approximate  model  of  the  NASA  gear  test  facility  in  Figure  5.2; 
for  gearbox  details,  refer  to  Chapter  IV.  The  intent  of  the  SEA  analysis  is  to  predict 
asymptotic  casing  plate  response  and  radiated  sound  pressure  from  the  gearbox,  on  an 
order  of  magnitude  basis.  Accordingly,  the  following  assumptions  are  made: 

1 .  Consider  the  case  of  4  identical  rolling  element  bearings  and  2  identical  shafts 
carrying  gears. 

2 .  The  vibratory  source  associated  with  the  static  transmission  error  excitation  is  at  the 
gear  mesh.  The  net  exchange  of  energy  through  the  gear  mesh  coupling  is 
neglected  and  the  vibratory  energy  is  assumed  to  flow  outward  from  the  gear  mesh 
source  to  the  casing  through  shafts  and  bearings.  The  source  is  assumed  to  be  far 
away  from  the  bearings. 

3.  Only  shaft  bending  vibration  is  coupled  to  the  casing  plate  flexural  motion. 
Torsional  modes  are  not  included  since  they  are  relatively  low  in  number  compared 
to  the  shaft  bending  modes.  In  steady-state,  both  shafts  have  equal  amount  of 
vibratory  energy  (Es)  associated  with  their  bending  motions  which  are  uncoupled 
from  the  external  motor  and  load. 

4.  The  portion  of  torsional  vibratory  energy  which  is  transmitted  to  the  load  is 
incorporated  in  the  dissipation  loss  factor  (rjs)  of  the  shaft-bearing  subsystem. 
Here  rjs  is  assumed  to  be  given  by  T|s=Ync^Tls(a>)  where  T|c*rjc(a>)  is  the  loss 


Mxsa(°J’z=°  )  =  Mxsa 

vysa(w,z=0")=0 


vysa(0)’  =Ls)  =  0 

^(£0,z=Lt  )=0 


Bearing 


Flexible  Shaft 


Bearing 


Figure  5. 1 6  E^froplccase  IV:  (a)  a  generic  geared  rotor  system  with  casing  and  mounts 
and  (b)  boundaiy  conditions  for  a  finite  shaft-bearing  system. 
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factor  of  the  casing;  y  is  expected  to  be  higher  than  one  due  to  additional  losses  at 
the  bearings  and  load. 

5 .  Bearing  dimensions  are  assumed  to  be  much  smaller  than  the  wavelength  on  the 
casing  plate.  Also,  bearings  are  not  expected  to  be  exchanging  vibratoiy  energy. 

6.  Modal  densities  of  shaft  tig  and  casing  plate  n^.  are  given  by  equation  (5.1). 

7.  A  diffuse  vibration  field  is  assumed  for  the  external  casing-mount  system.  The 
energy  dissipation  mechanism  for  this  system  includes  acoustic  energy  radiated 
from  the  plate  in  addition  to  the  structural  damping. 

Application  of  the  SEA  principle  will  result  in  2  subsystems,  one  internal 
consisting  of  2  shafts  and  4  bearings  and  one  external  consisting  of  2  casing  plates  and 
mounts  similar  to  Figure  5.10b.  The  internal  subsystem  in  this  case  can  store  the  total 
vibratory  energy  of  2  shafts.  Vibratory  energy  transfer  as  viewed  through  Figure  5.10b 
represents  the  algebraic  sum  of  the  energy  transfers  through  4  bearings. 


5.7.2  Coupling  Loss  Factor 

The  coupling  loss  factor  qsc  between  the  internal  (shafts-bearings)  and  external 
(casing)  systems  is  derived,  based  on  the  formulation  given  by  equations  (5.11)  and 
(5.12).  This  implies  that  vibratory  energy  transfer  is  associated  with  only  the  dynamic 
moment  at  the  bearing,  and  the  casing  plate  is  infinite  with  respect  to  the  bearing 
dimensions.  The  driving  point  impedance  of  the  shaft  Zs  is  derived  next  given  the 
boundary  conditions  shown  in  Figure  5.16b.  At  the  left  bearing  end  (z=0*), 
vyga(<o,z=0')=0,  and  Mxsa(co,z=0' )=Mxsa  which  are  equivalent  to  Mxsa(co,z=0+)=Mxsa, 
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and  FySa((O)O+)-(kbyexA5ex0x)Mxsa+vysa(co(O+)  { ^byy  }  /0®)* 
The  other  bearing  end  at  z=L*  with  zero  translational  and  angular  velocities  is  described 
at  z=L’  as 

(kb0x0x  (<°>  L  s  )  +^by 0 xvysa( ^*^5) ) /0M)  » 

Fysa(co,z=Ls)=— (kby0x  ^  (to,Ls)  +  kbyyvysa(®»^g))/(*®)  (5.21) 


Using  equations  (5.5)  and  (5.13),  the  above  boundary  conditions  are  evaluated  to 
formulate  an  algebraic  problem  of  the  type  given  by  equation  (5.14).  The  nonzero 
elements  of  the  coefficient  matrix  [B]  are 
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B42  =  <-ik!tby8>-|£byy  +  iElk*>e't'L‘  ! 

B43  =  «[skbyex-kbyy+EIk*)rk'1''  ; 

B44  =  (-k4kbye>-kbW-EIk»)'k'L'  <5'22> 

The  right  hand  side  vector  { b }  of  equation  (5.14)  is  Mysa  {  kby0x/kb0xex, 
2 

io)/(EIk^),0,0  )  ,  and  the  shaft  impedance  Zs  is  evaluated  numerically  using  equation 
(5.8). 


5.7.3  Vibroacoustic  Response 

Consider  the  total  energy  dissipation  by  the  casing  as  the  sum  of  the  energy  terms 
associated  with  structural  damping  dissipation  and  energy  escape  via  sound  radiation. 
Hence  the  total  dissipation  loss  factor  T|cj  for  the  casing  of  area  A,,  and  mass  n^  is 


t\cT«»)  =  i\c  + 


(omc 


(5.23) 


where  zG  is  the  characteristic  impedance  of  the  surrounding  medium,  and  CTC  is  the 
radiation  efficiency  of  the  casing.  In  this  analysis,  two  radiation  efficiency  models  are 
used:  (i)  oci  =  1.0  for  an  ideal  radiator  as  for  many  gearboxes  the  measured  radiation 
efficiency  has  been  found  to  be  close  to  unity  [78],  and  (ii)  oc2  for  a  simply-supported 
rectangular  plate  [79].  Several  investigators  have  used  oC2  successfully  in  SEA 
applications  [60,66,67]. 
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Since  this  SEA  model  is  very  similar  to  the  one  developed  for  example  case  III, 

solutions  given  by  equations  (5.19)  and  (5.20)  are  still  valid  here  provided  the  existence 

of  multiple  paths,  as  opposed  to  only  one  path  in  example  case  IE,  is  recognized.  The 
mean-square  velocity  levels  of  the  shafts-bearings  system  <Vs>  and  casing-mount 
,r2 

system  <VC>  are 


<  V,  >  (<a)  = 


<V^>(o»  = 


ns(TlcT+4tlcs) 


ms®(TlsncT+4TlcTV+4llsTlcs)  ’ 

_ 4nsnsc _ 

mc“(TisncT+4ncTnsc+4TisTics) 


(5.24) 


Power  injected  into  the  system  is  developed  by  examining  the  internal  static  transmission 
error  excitation  e(t)=eh(<0)eia,t-  F°r  the  spur  gear  pair,  the  gear  mesh  frequency 
G)(rad/sec)=NgC2zM7C/30  where  Ng  is  the  number  of  gear  teeth  and  QzM(rpm)  is  the 
driven  shaft  speed.  By  definition,  the  input  power  by  the  gear  mesh  elastic  force 
Fh(0=khe(t)=eh(w)ei<nt  is  given  by  ris(co)=l/2  (k^)2  Re  { (l/Z^co))* } .  At  a  very  high 
frequency,  it  is  reasonable  to  assume  that  the  boundaries  of  the  shafts  do  not  affect  the 
gear  mesh  source  regime.  Hence  the  infinite  shaft  assumption  should  hold. 
Accordingly,  Zg  can  be  obtained  analytically  using  the  driving  point  shaft  impedance 
given  in  Reference  [77];  however,  the  gear  mass  mg  must  be  included 

Zs(co)  =2psAscs  (1 +  i)  +  itomg  (5.25) 

Sound  power  W  radiated  from  the  casing  is  computed  using  the  following  equation 
where  Ac  is  the  total  casing  surface  area 
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W(co)  =  z0A  c<  Vj(co)  >  o  c(co)  (5.26) 

Assuming  a  source  directivity  Q(co)  associated  with  geared  system  mounting  condition, 
the  sound  pressure  level  Lp  (dB  re  20|iPa)  in  the  far  field  at  distance  r  from  the  casing  is 

0(0))  4 

L  p(o»  =  L  W(C0)  +  10  log  10  (^2  +  W  <5-27> 

where  R(ce>)=aS/(l-a)  is  the  room  constant,  a  (to)  is  the  average  absorption 
coefficient,  S  is  the  room  surface  area  and  L\y  (dB  re  l.OpW)  is  the  sound  power  level. 

5.7.4  Experimental  Validation  and  Parametric  Studies 

The  SEA  formulation  for  this  example  case  is  verified  by  comparing  results  with 
vibroacoustic  responses  measured  on  the  NASA  gearbox  [55].  A  detail  description  of 
the  experimental  setup  and  its  system  parameters  have  been  given  in  Chapter  IV. 
Additionally  unweighted  Lp  at  r  «  0.38  m,  directly  above  the  surface  of  the  top  cover 
casing  plate  as  shown  in  Figure  5.2,  has  been  measured  with  a  B&K  type  2230  sound 
level  meter.  Predicted  and  measured  Lp  are  compared  in  Figure  5.17a  for  y=10,  oc=oC2 
and  bearing  preload  5zm=0.04mm.  The  second  term  in  equation  (5.27)  is  dominated  by 
Q/(4itr2)  since  4/R  «  Q/^itr2)  in  this  case  due  to  a  «  1 .0  and  room  surface  area  S 
being  very  large.  Predictions  are  found  to  be  within  ±10dB  of  the  measured  values  for 

typical  structural  dissipation  loss  factor  0.004£r|c£0.04.  Figure  5.17b  shows  the 
comparison  of  spatially  averaged  mean-square  casing  acceleration  level  (<o2<V2>) 

between  theory  and  experiment.  Experimental  curve  represents  the  averaged  value  of 
the  measurements  made  at  3  casing  plate  locations  (top  plate,  side  plate  with 


Casing  Acceleration  (g)  Sound  Pressure  Level  (dB) 
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1/3  Octave  Band  Center  Frequency  (Hz) 

(a) 


1/3  Octave  Band  Center  Freqeuncy  (Hz) 


(b) 

Figure  5.17  Comparison  between  theory  and  NASA  experiment  [55]  for  example  case 
IV  (y=10,  Q=2,  Cc=ac2, 6zrn=0.04mm).  (a)  Sound  pressure  level  Lp. 
(b)  Spatially  averaged  casing  acceleration. 
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bearings  and  side  plate  without  bearings).  Again,  predicted  acceleration  spectra,  which 
are  similar  to  that  found  for  the  sound  level,  are  found  to  be  in  reasonable  agreement 
with  the  measured  spectrum  given  the  appropriate  values  of  T|c  and  y.  In  general,  we 
observe  that  the  response  level  decreases  with  increasing  T|c.  Also,  comparisons 
suggest  that  rjc  in  this  system  may  be  frequency  dependent. 

Next  we  vary  y  but  keep  qc=0.02  constant  in  Figure  5.18.  Comparison  between 
theory  and  measurements  also  indicates  that  y=10  is  the  best  fit  for  the  experimental  data 
especially  at  the  higher  frequencies.  Here,  radiation  efficiency  0cj  and  bearing  axial 
preload  5zfn=0.04mm  have  been  specified. 

Now  we  investigate  the  effects  of  casing  plate  radiation  efficiency  oc  and  bearing 
preload  or  mean  axial  displacement  8zm  on  Lp.  Figure  5.19  compares  acl  and  cc2. 
Based  on  the  comparison  with  experiment  it  seems  that  the  radiation  efficiency  of  the 
NASA  gearbox  is  better  modeled  with  CTcl  although  the  measured  Lp  curve  is  mostly 
between  crci  and  <rc2  curves.  Differences  between  acl  and  <Jc2  are  significant  at  lower 
frequencies,  but  the  variation  never  exceeds  lOdB  over  the  entire  frequency  range  of 
interest.  It  may  be  noted  that  since  the  acoustic  energy  radiated  W(to)  is  significantly 
smaller  than  the  energy  dissipated  by  the  system,  virtually  no  change  is  found  in  the 
predicted  casing  acceleration  spectra  by  varying  oc.  The  effect  of  bearing  preload  on  Lp 
is  shown  in  Figures  5.20a  and  5.20b  for  cC2  and  oci  respectively.  Figure  5.20 
indicates  that  Lp  is  lower  when  the  bearing  preload  is  increased.  A  similar  effect  is  seen 
for  the  casing  acceleration  level  which  is  consistent  with  the  deterministic  model 
prediction  of  Chapter  IV. 
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Figure  5.19  Sound  pressure  level  prediction  with  two  <TC  formulations  for  example 
case  IV:  ocj  =  1.0  (ideal  radiator)  and  aC2  for  simply-supported  plate 
(Y=10,  Q=2,  t)c=0.02,  Szm=0.04mm). 
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5.8  CONCLUDING  REMARKS 

The  vibration  transmission  through  bearing  has  been  analyzed  using  the  SEA 
technique.  A  new  procedure  has  been  developed  to  compute  the  coupling  loss  factor 
which  relies  on  the  solution  of  the  boundary  value  problem  at  the  bearing-casing 
interface.  This  scheme  incorporates  the  bearing  stiffness  matrix  developed  earlier  as  a 
part  of  the  deterministics  vibration  models  in  Chapters  II-IV.  Experimental  validations 
verify  the  proposed  theory  even  though  a  very  simple  model  for  the  geared  rotor  system 
has  been  presented.  A  more  detailed  SEA  model  including  energy  sharing  between 
subsystems  is  required  to  analyze  this  system.  This  is  left  for  future  research. 


CHAPTER  VI 


CONCLUSIONS  AND  RECOMMENDATIONS 

6.1  SUMMARY 

A  new  precision  rolling  element  bearing  model  is  proposed  for  the  analyses  of  the 
vibration  transmission  through  bearings  and  overall  system  dynamics  in  rotating 
mechanical  equipment  including  geared  rotor  systems.  Current  bearing  models,  based 
on  either  ideal  boundary  conditions  or  purely  translational  stiffness  element  description, 
can  not  explain  how  the  vibratory  motion  may  be  transmitted  from  the  rotating  shaft  to 
the  casing  and  other  connecting  structures  in  rotating  mechanical  equipment.  For 
example,  a  vibration  model  of  a  rotating  system  based  upon  the  existing  bearing  models 
can  only  predict  purely  in-plane  type  motion  on  the  flexible  casing  plate  given  only  the 
bending  motion  on  the  shaft.  However,  experimental  results  have  shown  that  the  casing 
plate  motion  is  primarily  flexural  or  out-of-plane  type.  Chapter  II  clarifies  this  issue 
qualitatively  and  quantitatively  by  developing  a  new  mathematical  model  for  the 
precision  rolling  element  bearings  from  basic  principles.  A  comprehensive  bearing 
stiffness  matrix  [K]t,m  of  dimension  6  is  proposed  which  clearly  demonstrates  a 
coupling  between  the  shaft  bending  motion  and  the  flexural  motion  on  the  casing  plate. 
A  numerical  scheme  which  involves  a  solution  to  nonlinear  algebraic  equations  is 
proposed  for  the  estimation  of  the  stiffness  coefficients  given  the  mean  bearing  load 
vector.  And,  a  second  method  which  requires  the  direct  evaluation  of  these  stiffness 
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coefficients  given  the  mean  bearing  displacement  vector  is  also  discussed.  Some  of  the 
translational  stiffness  coefficients  of  the  proposed  bearing  matrix  have  been  verified 
using  available  analytical  and  experimental  data.  Further  validation  of  [K]bm  is  not 
possible  as  coupling  coefficients  are  never  measured.  Also,  parametric  studies  on  the 
effect  of  unloaded  contact  angle,  preload,  or  bearing  type  are  included.  These  results 
lead  to  a  complete  characterization  of  the  bearing  stiffness  matrix. 

Chapter  III  extends  the  proposed  bearing  stiffness  formulation  and  demonstrates 
its  superiority  over  existing  models  in  vibration  transmission  analyses  for  a  generic 
single  shaft-bearing-plate-mount  system.  The  bearing  stiffness  matrix  [K]bm  is 
incorporated  in  discrete  system  models  using  lumped  parameter  and  finite  element 
modeling  techniques.  Shaft,  plate  and  mount  flexibilities  are  also  included  in  such 
models.  The  stability  issue  associated  with  the  proposed  bearing  model  is  addressed 
analytically  using  Liapunov's  stability  method  Mid  the  system  is  found  to  be  dynamically 
stable  provided  the  preloads  are  sufficiently  high.  Eigensolution  and  forced  harmonic 
response  to  the  following  rolling  element  bearing  system  example  cases  are  obtained 
using  our  formulation  and  results  are  compared  with  the  predictions  yielded  by  the 
current  vibration  models:  I.  rigid  shaft  and  plate  system  freely  suspended,  II.  rigid  shaft 
and  plate  supported  on  flexible  mounts,  and  IE.  an  experimental  setup  consisting  of  a 
flexible  shaft,  two  ball  bearings,  a  rectangular  plate  and  the  supporting  structure. 
Analytical  results  indicate  that  our  proposed  model  is  indeed  capable  of  predicting  plate 
rigid  body  angular  motion  or  plate  flexural  motion  as  excited  by  shaft  motion.  Such 
predictions  are  not  observed  in  simple  vibration  models.  Also,  lower  degrees  of 
freedom  models,  developed  by  several  previous  investigators  tend  to  underestimate  the 
resonant  frequencies  and  force  or  moment  transmissibilities  as  compared  with  our  multi- 
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degree  of  freedom  models.  Comparisons  between  our  model  and  the  available 
experiments  have  been  found  to  be  reasonably  good. 

The  overall  dynamics  of  a  geared  rotor  system  which  includes  a  spur  gear  pair, 
shafts,  rolling  element  bearings,  prime  mover  and  load  (attached  to  the  geared  rotor 
system  through  flexible  torsional  couplings),  rigid  or  flexible  casing,  and  compliant  or 
massive  mounts  is  studied  in  Chapter  IV.  Linear  time-invariant,  discrete  dynamic 
models  of  a  generic  geared  rotor  system  with  proportional  viscous  damping  are 
developed  using  lumped  parameter  and  dynamic  finite  element  techniques  which  are  then 
used  to  predict  the  vibration  transmissibility  through  bearings  and  mounts,  casing 
vibration  motion,  and  dynamic  response  of  the  internal  rotating  system.  The  proposed 
bearing  formulation  is  also  included.  Each  rotating  shaft  is  modeled  as  an  Euler  beam  in 
the  lumped  parameter  model  and  as  a  Timoshenko  beam  in  the  dynamic  finite  element 
model,  but  the  gyroscopic  moment  is  not  included.  Eigensolution  and  forced  harmonic 
response  studies  due  to  rotating  mass  unbalance  or  kinematic  transmission  error 
excitation  for  the  following  example  cases  are  obtained  using  our  formulation  and 
compared  with  simple  models  currently  available  in  the  literature  and/or  experiment:  I.  a 
single-stage  rotor  system  with  flexibly  mounted  rigid  casing  consisting  of  two  bearings 
as  a  special  case  of  the  geared  rotor  system,  n.  a  spur  gear  pair  drive  supported  by  four 
bearings  installed  in  a  flexibly  mounted  rigid  casing,  and  III.  an  experimental  setup 
consisting  of  high  precision  gear  and  pinion,  and  four  identical  rolling  element  bearings 
contained  in  a  flexible  casing  mounted  rigidly  on  a  massive  foundation.  In  example 
cases  II  and  01  of  Chapter  IV,  the  gear  mesh  stiffness  is  assumed  to  be  linear  and  time- 
invariant  Analytical  predictions  show  that  our  theory  is  indeed  capable  of  predicting 
bearing  and  mount  moment  transmissibilities  in  addition  to  the  force  transmissibilities. 
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Also,  flexural  vibrations  of  the  casing  plate  are  predicted  well  as  our  theory  is  in  good 
agreement  with  measurements  made  on  case  in  of  Chapter  IV;  such  predictions  are  not 
seen  in  simple  models. 

In  Chapters  II-IV,  the  deterministic  vibration  models  of  geared  rotor  system  with 
proposed  bearing  stiffness  matrix  incorporated  have  been  shown  to  be  reliable  up  to 
moderately  high  frequencies.  But  these  models  are  inadequate  at  very  high  frequencies 
when  the  modal  density  is  high.  Classical  vibration  models  usually  do  not  predict 
modes  accurately  in  this  frequency  regime,  require  large  computational  effort  and 
produce  spectra  at  many  spatial  points  which  are  difficult  to  interpret.  Accordingly,  we 
have  used  the  statistical  energy  analysis  (SEA)  method  to  predict  the  mean-square 
vibratory  response  of  internal  and  external  subsystems  in  a  geared  rotor  system.  The 
feasibility  of  applying  SEA  to  this  system  is  also  investigated  by  performing  modal 
analysis  of  a  gearbox.  Four  example  cases  are  analyzed  using  SEA;  I.  a  plate- 
cantilevered  rectangular  beam,  II.  case  I  with  circular  shaft-bearing  system  replacing  the 
cantilevered  rectangular  beam,  HI.  a  circular  shaft-bearing-plate-mount  system,  and  IV  a 
simple  geared  rotor  system.  In  the  first  two  example  cases,  we  have  revised  and 
extended  Lyon  and  Eichler's  plate-cantilevered  rectangular  beam  problem  [56,58]  to 
improve  the  coupling  loss  factor  prediction  and  to  formulate  the  vibratory  energy 
transfer  problem  through  rolling  element  bearings.  The  third  and  final  example  cases 
compute  the  system  response  spectra  and  compare  them  with  measurements.  Good 
agreement  is  found  between  theory  and  experiment  provided  proper  values  of  the 
dissipation  loss  factors  and  bearing  preloads  are  used.  The  NASA  gearbox  radiation 
efficiency  is  found  to  be  nearly  unity  at  higher  frequencies. 


62  FUTURE  RESEARCH 

Several  areas  of  potential  research  problems  based  on  the  present  study  of  the 

vibration  transmission  through  bearings,  are  identified  as  follows: 

1.  Analyze  vibration  transmission  through  hydrodynamic  bearings,  using  the 
proposed  work  on  rolling  element  bearings.  Develop  an  experimental  methodology 
to  estimate  the  bearing  stiffness  matrix  and  transfer  properties.  Also,  develop 
bearing  diagnostic  techniques  using  vibration  transmission  theory  proposed  in  this 
dissertation. 

2 .  Extend  the  proposed  overall  geared  rotor  system  vibration  model  with  spur  gears  to 
helical,  bevel  and  worm  gear  drives.  Develop  vibration  models  using  the  proposed 
bearing  model  to  analyze  multiple  transmission  paths.  Improve  gear  diagnostic 
procedures  using  analytical  bearing  transfer  properties. 

3 .  Refine  the  current  statistical  energy  analysis  formulation  of  a  generic  geared  rotor 
system  to  investigate  the  vibratory  energy  transfer  among  smaller  subsystems  and 
to  include  the  effects  of  rotating  shaft  torsional  modes.  Generate  analytical  and 
experimental  schemes  to  predict  coupling  loss  factors  for  typical  shaft-bearing-plate 
system.  Finally,  establish  gearbox  design  methodology  for  reduced  airborne  and 
structure-borne  noise. 
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